楊先鳳，李映潔，賴俊良，彭 博

（西南石油大學計算機科學學院，成都610500）

1 概述

超聲彈性成像是一種基于傳統(tǒng)超聲成像系統(tǒng)的非侵害性影像技術，它能夠將局部組織應變映射成彩色編碼的圖像信息，進而幫助臨床醫(yī)生確定組織硬度及相應的病理情況［1］。但是彈性成像系統(tǒng)在處理時涉及大量復雜運算，使其難以在臨床實時系統(tǒng)中得到廣泛應用［2］。

為了提高超聲彈性成像的質量，同時實現實時成像，本文提出一種基于GPU并行粒子群優(yōu)化的超聲彈性實時成像算法，并通過實驗分析該算法。

2 背景介紹

近年來，具有實時彈性成像能力的算法與方法得到了超聲工業(yè)界更多的關注與應用，如相位零估計（Phase Zero Estimation，PZE）［3］、聯合自相關［4］、聯合塊匹配與光流方法［5］、動態(tài)規(guī)劃［6］和基于先驗估計的互相關［7］等方法都已經部署在了商業(yè)化的超聲成像系統(tǒng)上。這些方法雖然得到了實時幀率，但是它們在提高計算速度的同時又導致了彈性圖質量降低等副作用。在較理想條件下，歸一化互相關（Normalized Cross Correlation，NCC）算法通?？梢蕴峁└哔|量彈性圖［8］，但是在臨床條件下，位移估計點易受噪聲干擾，因此這種采用孤立估計點計算位移的方式無法滿足更高精度的位移估計和生成高精度彈性圖的要求。粒子群優(yōu)化（Particle Swarm Optimization，PSO）算法是一種基于群智能的演化計算技術［9］，該算法通過個體之間的共享機制，使整個集群在解空間中從無序到有序的運動過程，搜索出復雜空間中的最優(yōu)解。PSO的優(yōu)點是流程簡單、易于實現，算法參數簡單，不需要繁雜的調整，同時PSO算法本身具有并行計算的特點。由于大量實驗已經證實，對求解最優(yōu)化問題，粒子群優(yōu)化算法具有收斂效率高、計算速度快的特點［10－11］。因此，近年來PSO算法被廣泛地應用于科學計算、工程應用和社會經濟等領域。對于超聲彈性的運動追蹤和位移估計過程而言，也可看作是一個最優(yōu)化問題。其目的是尋找能夠對壓縮前后的RF信號的運動情況進行準確描述的位移場。文獻［12］提出一個基于動態(tài)規(guī)劃算法的超聲彈性成像方法，該方法設計了一個結合回波幅度相似性與位移連續(xù)性的代價函數，并通過動態(tài)規(guī)劃算法優(yōu)化得到位移場。然而這個方法得到最終的位移場需要多層動態(tài)規(guī)劃方法或解析最小化方法［6］才能獲得。本文提出一個基于粒子群優(yōu)化的超聲彈性成像方法，與動態(tài)規(guī)劃多層優(yōu)化不同，基于粒子群優(yōu)化算法可以直接得到精確的位移。

隨著圖形處理器（Graphics Processing Unit，GPU）技術的快速發(fā)展，當前GPU已經具有很強的并行計算能力，浮點計算能力可以達到同代CPU的10倍以上［10］。GPU中大部分資源都被設計用來進行數據處理，只有極少數用于數據緩存和指令流的控制［13］。所以在信號和圖像處理方面，相對于多核CPU，帶有多處理器的GPU具有更為明顯的優(yōu)勢［14］。統(tǒng)一計算設備架構（Compute Unified Device Architecture，CUDA）［15］是由 NVIDIA 公司提出的一種基于GPU的并行計算架構，它的推出可以更加方便地在GPU上實現各種領域的科學計算。在CUDA平臺下對本文方法進行GPU并行框架設計，使其最終能達到實時的超聲彈性成像效果。

3 基于粒子群優(yōu)化算法的超聲彈性成像

3．1 粒子群優(yōu)化算法

粒子群優(yōu)化算法利用搜索空間的位置代表候選解，將搜索和優(yōu)化過程模擬成粒子的位置移動更新過程，每個粒子在搜索時都具有初始速度和位置。在搜索過程中，不僅要考慮自己搜索的歷史最佳位置，又要關注領域內其他粒子的歷史最佳位置，在此基礎上進行位置的變化［16］。將該求解問題的目標函數適應度值用來確定當前粒子的位置好壞，將個體的適者生存過程類比為搜索和優(yōu)化過程中用較好的可行解來取代較差的可行解的迭代過程。

假設在一個D維空間中，有N個粒子，記粒子i的位置為xi＝（xi1，xi2，…，xiD），其中，xid∈［l，u］，d∈［1，D］，l和u分別為D維空間的下界與上界。將xi代入適應度函數f（xi）求出適應度值；記粒子i的速度為vi＝（vi1，vi2，…，viD），用戶指定的Vmax限制粒子速度在一定范圍內。記粒子i個體經歷過的最好位置為pi＝（pi1，pi2，…，piD）。標準的粒子群算法采用如下公式對粒子速度和位置進行更新［17］：

其中，學習因子c1和c2是非負常數，使粒子具備自我歸納和向群體中良好個體學習的能力，從而向自己的歷史最優(yōu)點和群體內或領域內的歷史最優(yōu)點接近，一般取c1＝c2＝2；r1和r2是隨機數，服從［0，1］上的均勻分布；pgbest是粒子群所有粒子中的最好位置；t為迭代次數，一般將最大迭代次數或粒子群迄今為止搜索到的最優(yōu)滿足預先設定的最小閾值作為迭代中止條件。

3．2 超聲彈性成像運動追蹤與位移估計

假定從壓縮組織變形前后得到的RF信號為RF1和RF2。RF1和RF2的尺寸為m×n（m表示一條回波線上采樣點的數量，n表示回波線的數量）。這里，適應度函數由標準的回波幅度相似度計算方法構成：

當ER越小，說明2條回波線越相似。在實際的計算過程中，通常可采用一種類似于相關性計算的比較窗口來計算它們的相似性，并且為了避免搜索過程中陷入局部最優(yōu)，本文采用一種類似于金字塔模式的多層次窗口模式計算對應估計位置的相關性。通常，過高的窗口重疊率會產生蠕蟲噪聲［18］，采用一定的間隔來計算對應采樣點的位移可以滿足應用需求，同時也避免了過于密集的計算。因此，對于采樣點i對應的總代價函數如下式：

在基于粒子群優(yōu)化算法中，di是所求的待估計點（待估計窗口）位移，最優(yōu)di值是通過粒子群優(yōu)化獲得，因此RF2（k＋di，j）的取值必須通過插值方式取得。本文采用線性插值方法。整個總代價函數作為粒子群優(yōu)化算法中的適應度函數。

本文提出的基于粒子群優(yōu)化算法的超聲彈性成像位移估計的具體實現步驟如下：

步驟1將壓縮前后的2個超聲回波RF信號文件讀入到內存中。

步驟2每一條回波線上的采樣點個數為ROW，共有COL條掃描線，設定每一條回波線上的窗口大小為win個采樣點，窗口間隔為k，則可以得到待估計點的數量M。用粒子群優(yōu)化計算每一個待估計窗口的位移delaym，其中，m＝1，2，…，M。

步驟3利用傳統(tǒng)的互相關算法求得位于中心位置回波線的位移作為下一步粒子群優(yōu)化估計的引導位移，這可以減少粒子群優(yōu)化算法的搜索范圍。

步驟4隨機初始化每一個粒子的位置Xi和速度Vi。

步驟5從t＝1到最大迭代次數，執(zhí)行下面的循環(huán)：

（1）由式（3）計算出粒子的適應度值fi；

（2）計算粒子的個體最優(yōu)與群體最優(yōu)；

（3）通過粒子群優(yōu)化算法的速度Vi和位置Xi更新式（1）和式（2）對粒子狀態(tài)進行更新。

步驟6迭代結束以后，得到的群體最優(yōu)位置和引導位移的和即為待估計點的位移，回步驟4繼續(xù)下一個待估計點的位移估計。

3．3 并行超聲彈性成像框架

從粒子群優(yōu)化的適應度函數式（3）計算知道，對于每一條信號上每一個窗口的計算都是數據獨立的。因此，基于粒子群優(yōu)化算法的超聲彈性成像的并行性主要體現在以下2個方面：

（1）一條掃描線上的每個窗口的位移是可以并行計算的；

（2）每一條掃描線是相互獨立的，也是可以并行計算的。

因此，從并行算法的角度，每個待估計點的所在窗口都可以同時進行位移估計，也就是讓每一個GPU線程處理一個待估計點所在窗口的所有采樣點。進一步需要考慮時如何采用有效的線程結構和優(yōu)化的存儲器以提高帶寬的使用率。

在CUDA平臺上，warp是調度和執(zhí)行的基本單位，一個warp包含32個線程。為了防止過小的warp浪費計算資源，每個線程塊中線程的數量應該設置為warp尺寸的整數倍，即線程數量的設置至少需要滿足的條件是一個線程塊內的線程數是32的倍數，并且最好不低于64個。使用紋理存儲器是得到高帶寬利用率的手段之一，壓縮前后的2個I／Q數據幀存放在紋理存儲器中。紋理存儲器提供的緩存機制使對紋理緩存的拾取只需要一個存儲器時鐘，這樣可以達到較高的讀取數據?；诹Ｗ尤簝?yōu)化算法的超聲彈性成像并行實現的具體實施思路：在GPU中建立所需且有效的線程結構，并為每個待估計點開辟獨立的計算空間，同時保證合理的寄存器使用數量，每個線程運行中均會執(zhí)行PSO算法以針對特定待估計點進行位移尋優(yōu)。GPU的多核并行特性使多個窗口待估計點位移計算的時間縮短成與一個窗口計算的時間相似，從而提高運算效率。本文方法的GPU并行實現框架如圖1所示。

圖1 GPU并行框架

由于在仿真實驗中驗證了較小的種群規(guī)模與迭代次數能夠實現準確的超聲彈性成像的位移追蹤。因此，與傳統(tǒng)粒子群算法本身的并行不同［19－20］，本文將粒子群優(yōu)化算法實現為device函數，通過一個并行計算所有估計點的kernel調用即可滿足實時計算需求。基于CUDA的并行執(zhí)行kernel函數如圖2所示。

圖2 基于CUDA的kernel函數

其算法描述為：

其中，tx代表每一個線程，在＿＿global＿＿函數里面調用粒子群優(yōu)化算法函數計算每一個窗口的位移，由于使用＿＿global＿＿的函數，只可以在設備上執(zhí)行，因此要調用粒子群優(yōu)化函數PSO（），只能將PSO（）作為＿＿device＿＿函數才以在設備上調用，且只可以在設備上調用。PS0（）函數的具體步驟與上一節(jié)所述一樣。

4 實驗與結果分析

4．1 仿真實驗

仿真使用文獻［19］提出的模擬壓縮前后的生物組織變化數學模型。在本文中，仿真的散射子模型內部包含一個硬球，模型的寬度為3cm，深度為3cm，硬球的半徑為0．5cm，位于模型內部中間位置。這個模型由200 000個散射子構成，散射子的強度符合高斯分布。

該模型分別考慮到散射的橫向和縱向移動，模型中散射子的橫向和縱向移動位移定義如下：

其中，R是硬球的半徑；r是散射子到坐標中心的距離，并且r2＝x2＋z2；K是硬球的楊氏模量與背景的楊氏模量之比；P是外力應變值，接近等于背景的軸向應變；v（x，z）和u（x，z）分別是散射子的軸向位移和橫向位移，本文僅需要軸向位移信息，設定u（x，z）＝0即可。

4．1．1 種群規(guī)模與迭代次數對算法性能的影響

種群規(guī)模與迭代次數是影響粒子群優(yōu)化算法的優(yōu)化精度和收斂速度的2個重要參數，為了解決本文所提出的位移追蹤的優(yōu)化問題，在保證位移追蹤精度的同時，選擇適合的種群規(guī)模與迭代次數對提高算法的計算效率具有重要影響。

通過改變種群規(guī)模與迭代次數，評估粒子群優(yōu)化在尋找最佳適應度值（即最優(yōu)解）的性能優(yōu)劣，以及確定算法中所需要的迭代次數。

從圖3的實驗結果可以看出，在針對位移追蹤的優(yōu)化問題上，粒子群優(yōu)化算法表現出的特性為：隨著種群規(guī)模的增大，迭代次數變小，且算法優(yōu)化精度基本接近。所以PSO算法只需要較少的迭代次數就能獲得足夠的位移追蹤精度。

圖3 種群規(guī)模與迭代次數對算法性能的影響

4．1．2 位移估計準確性

為了評估粒子群優(yōu)化算法與傳統(tǒng)的互相關算法在位移估計上的偏差，本文分別采用互相關算法和粒子群算法計算位于中心位置的回波線的位移，通過均方根誤差方法對得到的位移與理想情況下的位移進行對比分析。

圖4顯示了基于粒子群優(yōu)化的運動追蹤和位移估計算法與傳統(tǒng)的互相關方法（拋物線插值）的位移估計精度。圖4（a）為整條回波線的位移比較，從中很難區(qū)分2種方法的優(yōu)劣，它們得到的位移估計值與真實位移值幾乎完全吻合。圖4（b）為對圖4（a）前100多個位移估計點放大圖，從圖中可以看出，粒子群優(yōu)化方法得到的位移與真實位移更為接近。顯示了粒子群優(yōu)化方法能獲得更好的位移估計精度。而互相關方法與粒子群算法所得到的整條線位移的均方根誤差分別為0．007 3和0．006 3。

圖4 算法真實位移的比較

4．1．3 計算速度

本文中CPU串行計算采用Intel Xeon E3－1220 CPU，GPU并行計算采用NVIDIA Tesla K10GPU卡實現，K10共有3 072個流處理器核和8GBDDR5內存。CUDA驅動版本為5．0，使用VS 2008作為編程工具。

表1為粒子群優(yōu)化算法的CPU串行實現與粒子群優(yōu)化算法的GPU并行實現的運行時間對比結果。為了公平起見，CPU和GPU實現都采用相同的參數設置，計算時窗口大小設置為20個采樣點，相鄰窗口的重疊部分為75%，種群數量為20，迭代次數為15次。表中加速比為本文所提方法的CPU串行實現與GPU并行實現的比較結果。表1中前3行的計算結果針對原始RF信號行不變，而回波線的數量成倍增加。而后3行的計算結果針對原始RF信號保持列不變，采樣點數量逐步成倍增加。從表1中可以得出，不管數據大小為多少，與本文方法的CPU實現相比，該GPU并行實現方法獲得了較好的加速比。針對行保持不變的數據，加速比有小幅度的變化；而針對列保持不變的數據，采樣點的數量對加速比卻有較大的影響。表1中基于GPU并行實現的計算時間和加速比驗證了本文方法的并行計算框架能夠通過GPU實現超聲彈性成像位移估計算法，在處理512×128的數據規(guī)模時其GPU實現處理時間僅為12 ms，可以滿足一般標準輸出條件下的實時超聲彈性成像的需要。

表1 CPU串行與GPU并行時間對比

4．1．4 仿真成像結果

圖5顯示了基于粒子群優(yōu)化的運動追蹤和位移估計算法與傳統(tǒng)的互相關方法（拋物線插值）在仿真數據情況下產生的應變圖。在圖5（a）中，得到的應變圖有幾處明顯的白色跳變區(qū)域。這說明互相關方法在計算位移時容易受到噪聲數據的干擾，這種干擾造成位移的錯誤計算，以至于應變估計放大了這種位移錯誤。圖5（b）顯示了粒子群優(yōu)化得到的應變圖，圖中無明顯跳變錯誤，這說明粒子群優(yōu)化具有較好的抗噪能力。

圖5 2種算法生成的應變圖

4．2 體模實驗

體模實驗采用專用于彈性成像研究的Model 049彈性體模出來的實時掃描射頻回波信號進行體模實驗。本實驗選擇的成像物體直徑為10mm、彈性模量為63kPa的硬包容物（Type IV）。實驗用探頭型號為SA5L38B的128陣元的線陣探頭，中心頻率為5MHz，75%的分數階帶寬。實驗使用的超聲系統(tǒng)是iMago C21超聲機，系統(tǒng)RF信號采樣頻率設為40 MHz。在數據采集過程中，未使用任何額外的控制設備情況下保持在一個恒定的速度將探頭進行徒手軸向壓縮／釋放。

圖6為對壓縮前后的2幀信號，分別采用互相關分析和粒子群優(yōu)化方法的串行和并行方法計算得到的彈性圖，從圖中可以得到，在真實體模數據情況下，粒子群優(yōu)化的2種方法都能獲得較好質量的彈性圖。

圖6 體模數據產生的彈性圖

圖6（c）所采用的并行粒子群優(yōu)化算法為了提高計算速度采用的計算精度是單精度，算法實現過程中的一些數學函數也采用的是低精度的版本，但是對于本文中的彈性成像的應用是足夠的。從視覺上來看，圖6（b）和圖6（c）2種方法幾乎沒有任何差別，這也說明文中所選擇的種群數量和迭代次數都是合適的。

5 結束語

本文將群智能算法中的粒子群優(yōu)化算法應用于超聲彈性成像中，通過CUDA平臺并利用GPU實現高效并行計算。實驗結果表明，本文算法不僅能夠獲得較好的位移質量，同時GPU并行實現的方法能夠保證彈性圖的實時處理。本文算法的提出對于拓展群智能算法的應用領域具有積極作用，同時對于超聲彈性成像中的最優(yōu)化問題求解也具有重要的現實意義，下一步將對本文算法繼續(xù)進行臨床數據測試，并優(yōu)化其GPU的并行實現性能。

［1］彭 博，諶 勇，劉東權．基于GPU的超聲彈性成像并行實現研究［J］．光電工程，2013，40（5）：97－105．

［2］張 霞，何興無．CUDA平臺下的超聲彈性成像并行處理算法［J］．計算機與數字工程，2012，（9）：113－116．

［3］Lindop J E，Treece G M，Gee A H，et al．3DElastography Using Freehand Ultrasound［J］．Ultrasound in Medicine ＆ Biology，2006，32（4）：529－545．

［4］Shiina T，Nitta N，Sjsum E U，et al．Real Time Tissue Elasticity Imaging Using the Combined Autocorrelation Method ［J］．Journal of Medical Ultrasonics，2002，29（3）：119－128．

［5］Zhou Yongjin，Zheng Yongping．A Motion Estimation Refinement Framework for Real－time Tissue Axial Strain Estimation with Freehand Ultrasound［J］．IEEE Tran－sactions on Ultrasonics，Ferroelectrics and Frequency Control，2010，57（9）：1943－1951．

［6］Rivaz H，Boctor E，Foroughi P，et al．Ultrasound Elastography：A Dynamic Programming Approach［J］．IEEE Transactionson Medical Imaging，2008，27（10）：1373－1377．

［7］Zahiri A R，Salcudean S E．Motion Estimation in Ultra－sound Images Using Time Domain Cross Correlation with Prior Estimates ［J］．IEEE Transactions on Bio－medical Engineering，2006，53（10）：1990－2000．

［8］Hoyt K，Forsberg F，Ophir J．Comparison of Shift Estimation Strategies in Spectral Elastography［J］．Ultra－sonics，2006，44（1）：99－108．

［9］Kennedy J，Kennedy J F，Eberhart R C．Swarm Intelligence［M］．［S．l．］：Morgan Kaufmann，2001．

［10］左顥睿，張啟衡，徐 勇，等．基于GPU的并行優(yōu)化技術［J］．計算機應用研究，2009，26（11）：4115－4118．

［11］王燕燕，葛洪偉，王娟娟，等．一種動態(tài)分組的粒子群優(yōu)化算法［J］．計算機工程，2015，41（1）：180－185．

［12］Rivaz H，Boctor E M，Choti M A，et al．Real－time Regularized Ultrasound Elastography［J］．IEEE Transactions on Medical Imaging，2011，30（4）：928－945．

［13］劉 丹，趙廣輝，夏紅霞，等．GPU加速分子動力學模擬的熱力學量提?。跩］．計算機應用研究，2010，27（5）：1820－1822．

［14］Zhao Jieyi，Tang Min，Tong Ruofeng．Connectivitybased Segmentation for GPU－accelerated Mesh Decompression［J］．Journal of Computer Science and Technology，2012，27（6）：1110－1118．

［15］都志輝，李三立．高性能計算之并行編程技術——MPI并行程序設計［J］．北京：清華大學出版社，2001．

［16］Spears W M，Green D T，Spears D F．Biases in Particle Swarm Optimization［J］．International Journal of Swarm Intelligence Research，2010，2（1）：34－57．

［17］張維存．蟻群粒子群混合優(yōu)化算法及應用［D］．天津：天津大學，2007．

［18］崔少國，劉東權．使用二維小波收縮法去除彈性成像蠕蟲噪聲［J］．生物醫(yī)學工程學，2011，8（3）：460－464．

［19］陳 風，田雨波，楊 敏．基于CUDA的并行粒子群優(yōu)化算法研究及實現［J］．計算機科學，2014，（9）：263－268．

［20］蔡 勇，李光耀，王 琥．基于CUDA的并行粒子群優(yōu)化算法的設計與實現［J］．計算機應用研究，2013，29（8）：2415－2418．