陳胤燃 何 瓊 羅建文*

1(清華大學(xué)醫(yī)學(xué)院生物醫(yī)學(xué)工程系，北京 100084)2(清華大學(xué)生物醫(yī)學(xué)影像研究中心，北京 100084)

基于GPU并行計(jì)算的超聲波束合成方法

陳胤燃1,2何 瓊1,2羅建文1,2*

1(清華大學(xué)醫(yī)學(xué)院生物醫(yī)學(xué)工程系，北京 100084)2(清華大學(xué)生物醫(yī)學(xué)影像研究中心，北京 100084)

逐列掃描式的聚焦成像方法會(huì)限制超聲成像幀頻的提高。采用平面波發(fā)射模式只需要一次發(fā)射和接收即可獲得一幅完整的超聲圖像，平面波成像是可以大幅度地提高成像幀頻，從而實(shí)現(xiàn)超聲超高速成像的方法之一。但是現(xiàn)有的超聲波束合成器無(wú)法達(dá)到超聲超高速成像對(duì)于計(jì)算能力的要求。對(duì)基于平面波成像的超聲延時(shí)-疊加波束合成算法進(jìn)行并行性分析，在此基礎(chǔ)上設(shè)計(jì)并實(shí)現(xiàn)兩種基于圖形處理器(GPU)并行計(jì)算的超聲平面波成像波束合成方法——基于2個(gè)Kernel和基于1個(gè)Kernel的并行波束合成方法。兩種方法的主要區(qū)別在于波束合成中對(duì)延時(shí)值的計(jì)算和存儲(chǔ)策略的不同處理，仿體實(shí)驗(yàn)證明兩種方法的計(jì)算幀頻分別達(dá)到2 178 和2 453 幀/s。相比于普通的方法，這兩種基于GPU的并行波束合成方法的計(jì)算幀頻分別提速99倍和111倍。實(shí)驗(yàn)結(jié)果表明，GPU波束合成器相比于傳統(tǒng)方法，可以大幅度提高計(jì)算能力。

超聲成像；波束合成；并行計(jì)算；圖形處理器

引言

傳統(tǒng)的超聲B模式成像利用多次發(fā)射和接收聚焦聲波來(lái)合成一幅完整的圖像，該成像模式下的理論最高幀頻一般在幾十到一百多幀每秒之間。如果采用多焦點(diǎn)發(fā)射的模式，該理論最高幀頻將相應(yīng)降低，降低的倍數(shù)等于多焦點(diǎn)增加的個(gè)數(shù)。為了進(jìn)一步提高超聲成像的幀頻，隨著超聲成像系統(tǒng)軟硬件水平的不斷提高，超高速超聲成像逐漸發(fā)展為如今超聲成像的研究熱點(diǎn)之一。

超高速超聲成像的概念最早由Bruneel等在1977年提出[1]。在1979年，和Bruneel同一個(gè)組的Delannoy等開(kāi)發(fā)出了一個(gè)超聲系統(tǒng)，可以通過(guò)一次的脈沖發(fā)射，利用模擬域的并行數(shù)據(jù)處理而獲得一幅圖像，該系統(tǒng)可以達(dá)到1 000幀/s的速度(70條線)[2]。在Bruneel和Delannoy等之后，1984年，Shattuck等在相控陣成像中通過(guò)在一次發(fā)射后同時(shí)計(jì)算多條B模式掃描線的方法來(lái)提高成像速度[3]。1991年，Lu和Greenleaf提出了基于非衍射波束的超聲成像來(lái)提高成像速度[4]。從20世紀(jì)90年代開(kāi)始，F(xiàn)ink等提出了平面波發(fā)射和并行波束合成的方法來(lái)實(shí)現(xiàn)超高速成像，該方法可以將幀頻提高到數(shù)千幀每秒[5-7]。平面波成像在發(fā)射過(guò)程中同時(shí)激勵(lì)超聲探頭的所有陣元，發(fā)出波前形狀為平面的脈沖超聲波，平面波在傳播過(guò)程中覆蓋整個(gè)成像區(qū)域。在接收過(guò)程中，各個(gè)陣元接收聲場(chǎng)中的回波信號(hào)，經(jīng)過(guò)并行波束合成和其他后處理之后得到一幅完整的B模式圖像，平面波成像只需要一次發(fā)射/接收過(guò)程即可以獲得一幅圖像，因此其理論最高成像幀頻為

(1)

式中，Z是成像區(qū)域的深度，c是超聲波在組織中傳播的速度。

通常認(rèn)為c的值為1 540m/s，假設(shè)圖像深度為5cm，那么對(duì)應(yīng)的理論最高幀頻為15 400f/s(幀/s)超高的幀頻可以帶來(lái)絕佳的時(shí)間分辨率，從而能夠捕捉到組織運(yùn)動(dòng)中的更多細(xì)節(jié)，有效地輔助疾病的診斷，如實(shí)現(xiàn)剪切波彈性成像、超高速超聲多普勒成像等。

就超聲波發(fā)射/接收的環(huán)節(jié)而言，平面波成像相比于傳統(tǒng)聚焦成像能夠大幅度地提高幀頻，但要真正實(shí)現(xiàn)超高速成像，還需要解決成像過(guò)程中因?yàn)閱挝粫r(shí)間內(nèi)數(shù)據(jù)量大幅度增加而帶來(lái)的計(jì)算能力不足的問(wèn)題?，F(xiàn)有臨床超聲設(shè)備的波束合成環(huán)節(jié)都通過(guò)數(shù)字信號(hào)處理器(DSP)、現(xiàn)場(chǎng)可編程門陣列(FPGA)等設(shè)備來(lái)實(shí)現(xiàn)，這些設(shè)備的計(jì)算能力還無(wú)法滿足超聲超高速成像波束合成實(shí)時(shí)計(jì)算的要求。

圖形處理器或圖像處理單元(graphicsprocessingunit,GPU)，由于具有強(qiáng)大的數(shù)值計(jì)算能力和高度的可并行性，使得它成為通用計(jì)算領(lǐng)域的強(qiáng)大工具，也是實(shí)現(xiàn)超聲超高速成像的工具之一。夏春蘭等利用GPU實(shí)現(xiàn)了從超聲射頻(RF)信號(hào)到B模式圖像的各項(xiàng)關(guān)鍵環(huán)節(jié)的處理，顯著地提高了成像的速度[8]。一些學(xué)者將GPU應(yīng)用于超聲的波束合成，并利用多個(gè)GPU實(shí)現(xiàn)了平面波成像波束合成的并行處理，在一定的成像參數(shù)下，實(shí)現(xiàn)了每秒數(shù)千幀的計(jì)算幀頻[9-12]。

本研究采用平面波成像方式和延時(shí)-疊加波束合成算法，利用一塊商用GPU，采用NVIDIA公司開(kāi)發(fā)的統(tǒng)一計(jì)算架構(gòu)(computeunifieddevicearchitecture,CUDA)，結(jié)合C語(yǔ)言進(jìn)行編程[13]，設(shè)計(jì)并實(shí)現(xiàn)了兩種基于GPU并行計(jì)算的平面波成像波束合成方法。

1 方法

1.1 CUDA編程模型

1.1.1 GPU并行計(jì)算模型

在介紹CUDA架構(gòu)的編程模型之前，首先了解幾個(gè)和CUDA編程有關(guān)的概念[13-16]：

1) Kernel：指可以在GPU中運(yùn)行的函數(shù)，也稱之為內(nèi)核函數(shù)，CUDA程序中需要交給GPU處理的命令都需要寫在Kernel函數(shù)中，供GPU調(diào)用。

2) Grid：是一個(gè)虛擬的概念，中文名稱“線程網(wǎng)格”或“網(wǎng)格”，它是一個(gè)和Kernel相對(duì)應(yīng)的概念，通常一個(gè)Kernel對(duì)應(yīng)著一個(gè)特定的Grid，即運(yùn)行一個(gè)Kernel都需要生成一個(gè)Grid。

3) Block：中文名稱為“線程塊”，Block是Grid的最小執(zhí)行單位，一個(gè)Grid中可以劃分出多個(gè)Block，各個(gè)Block在任務(wù)中是并行執(zhí)行的，Block之間沒(méi)有固定的執(zhí)行順序，也不能進(jìn)行通信，Grid里最多可以劃分出3個(gè)維度的Block。

4) Thread：中文名稱“線程”，是執(zhí)行計(jì)算任務(wù)的最小單位，一個(gè)Block里可以劃分出最多3個(gè)維度的Thread，同一個(gè)Block里的Thread在執(zhí)行任務(wù)的時(shí)候可以通過(guò)共享存儲(chǔ)器進(jìn)行通信。

圖1說(shuō)明了一個(gè)Kernel中的Grid、Block和Thread之間的抽象關(guān)系，可以看出在Kernel中存在著兩個(gè)層次的并行：Grid中Block之間的并行和Block中Thread的并行。

使用GPU進(jìn)行并行計(jì)算時(shí)，需要根據(jù)待處理數(shù)據(jù)的具體特征，合理設(shè)計(jì)Grid和Block的維度、Grid中每一個(gè)維度上的Block數(shù)目和Block中每一個(gè)維度上的Thread數(shù)目。

圖1 一個(gè)Kernel中Grid、Block和Thread之間的抽象關(guān)系(以二維為例)Fig.1 Abstract relations of grid, block and thread in one kernel (a 2-D case)

1.1.2 GPU數(shù)據(jù)存儲(chǔ)架構(gòu)

GPU中有多種類型的存儲(chǔ)器，它們和Grid、Block、Thread之間的讀寫訪問(wèn)關(guān)系如表1所示。

表1 GPU中各種存儲(chǔ)器的比較 [13-16]

利用GPU進(jìn)行并行計(jì)算時(shí)，數(shù)據(jù)在其中的存儲(chǔ)策略合理與否對(duì)計(jì)算的效率有很大的影響。在進(jìn)行波束合成并行計(jì)算時(shí)，需要盡可能地提高數(shù)據(jù)讀寫的速度。共享存儲(chǔ)器在所有的存儲(chǔ)器中，訪問(wèn)延遲較小，帶寬較高，可以同時(shí)被一個(gè)Block里的所有Thread并行訪問(wèn)。在并行波束合成的具體實(shí)現(xiàn)過(guò)程中，使用共享存儲(chǔ)器可優(yōu)化數(shù)據(jù)存儲(chǔ)策略，在一定程度上提高了波束合成的速度。

1.2 GPU并行波束合成方法設(shè)計(jì)和實(shí)現(xiàn)

1.2.1 延時(shí)-疊加波束合成算法并行性分析

采用延時(shí)-疊加波束合成算法，首先需要計(jì)算不同深度的延時(shí)值，根據(jù)延時(shí)-疊加算法的原理，圖像中同一深度上各個(gè)采樣點(diǎn)的延時(shí)模型是完全一致的。如圖2所示，第i行上的點(diǎn)L_i1和L_i2處在同一深度上，它們?cè)诟鱾€(gè)對(duì)應(yīng)通道中的延時(shí)完全一致(在圖中用相同的色塊表示)，可以共享2R+1個(gè)延時(shí)值，2R+1是波束合成接收孔徑的大小。因此，同一深度的各個(gè)采樣點(diǎn)可以通過(guò)共享一組延時(shí)值，實(shí)現(xiàn)并行合成。再分析同一通道不同深度的采樣點(diǎn)，可以發(fā)現(xiàn)不同深度的采樣點(diǎn)在進(jìn)行合成時(shí)是相互獨(dú)立的。同樣如圖2所示，第i行和第j行的采樣點(diǎn)在疊加計(jì)算時(shí)互不影響，也可以實(shí)現(xiàn)并行合成。不同深度點(diǎn)的延時(shí)值不同，在圖中由不同灰度的方塊表示。

圖2 延時(shí)-疊加波束合成算法并行性分析Fig.2 Analysis of parallel Computation in delay-and-sum beamforming

通過(guò)分析可知，延時(shí)-疊加波束合成算法至少存在著兩層可并行計(jì)算的環(huán)節(jié)。

1.2.2 2個(gè)Kernel實(shí)現(xiàn)并行波束合成

由于延時(shí)的計(jì)算和采樣點(diǎn)本身的值沒(méi)有關(guān)系，因此在波束合成前，可以先并行地計(jì)算出所有的延時(shí)，得到延時(shí)表，它的每一行代表了對(duì)應(yīng)深度采樣點(diǎn)所共享的一組延時(shí)值。在波束合成過(guò)程中，對(duì)于每個(gè)采樣點(diǎn)的所有延時(shí)不需要重復(fù)計(jì)算，而是通過(guò)“查找延時(shí)表”的方式，找到與之同一深度的所有延時(shí)值即可。

通過(guò)分析可以發(fā)現(xiàn)，延時(shí)表的計(jì)算過(guò)程存在著兩層的并行，因此可以將延時(shí)表的計(jì)算作為并行任務(wù)1。將任務(wù)1交給一個(gè)Kernel函數(shù)進(jìn)行計(jì)算，此處將該函數(shù)命名為Kernel1，通過(guò)計(jì)算得到的延時(shí)表存放于GPU中，供后面的疊加環(huán)節(jié)使用。

將查延時(shí)表進(jìn)行波束合成的環(huán)節(jié)作為并行任務(wù)2，并建立一個(gè)新的Kernel函數(shù)執(zhí)行，該函數(shù)此處命名為Kernel2，通過(guò)它讀取延時(shí)表的對(duì)應(yīng)數(shù)值進(jìn)行波束合成，Kernel1和Kernel2各自的Grid和Block設(shè)計(jì)及并行計(jì)算任務(wù)的劃分示意如圖3、4所示。

圖3 Kernel1對(duì)應(yīng)的Grid1和Block1設(shè)計(jì)及并行任務(wù)劃分Fig.3 Diagrammatic sketch of Grid1 and Block1 designs and task assignments in Kernel1

圖4 Kernel2對(duì)應(yīng)的Grid2和Block2設(shè)計(jì)及并行任務(wù)劃分Fig.4 Diagrammatic sketch of Grid2 and Block2 designs and task assignments in Kernel2

從圖3中可以總結(jié)出計(jì)算延時(shí)表的Kernel1過(guò)程中的兩層并行：

1) 延時(shí)表中第i行的2R+1個(gè)延時(shí)值，全部交給Grid1里的Block(0，i-1)進(jìn)行計(jì)算；

2) 該行的第k個(gè)延時(shí)值(共2R+1個(gè))，則交給Block(0,i)里的Thread(0,k-1)線程進(jìn)行計(jì)算。

從圖4可以總結(jié)出進(jìn)行疊加計(jì)算的Kernel2過(guò)程中的3層并行：

1) Grid2中第s行的M+1個(gè)Block負(fù)責(zé)計(jì)算圖像數(shù)據(jù)包中的第s張未合成圖像；

2) 這M+1個(gè)Block與第s張圖像的M+1行數(shù)據(jù)一一對(duì)應(yīng)，即Block(s, i)負(fù)責(zé)計(jì)算第s張圖像的第i行上的所有N+1個(gè)像素點(diǎn)的值；

3) 第s張圖像的第i行上的N+1個(gè)數(shù)據(jù)與Block(s,i)里的N+1條線程一一對(duì)應(yīng)，即Thread(0,j)負(fù)責(zé)合成第s張圖像的像素點(diǎn)P(i,j)。

1.2.3 1個(gè)Kernel實(shí)現(xiàn)并行波束合成

使用2個(gè)Kernel實(shí)現(xiàn)延時(shí)-疊加波束合成將 “計(jì)算延時(shí)”變成更加方便的“查找延時(shí)”，但是“查找延時(shí)”的設(shè)計(jì)存在著局限性。首先，查延時(shí)表的操作要求對(duì)應(yīng)于每一個(gè)接收孔徑的延時(shí)模式是完全一致的，這僅適用于平面波成像中的無(wú)角度偏轉(zhuǎn)發(fā)射/接收模式，若采用帶角度偏轉(zhuǎn)發(fā)射/接收的平面波成像方式，其延時(shí)模式更加復(fù)雜，因此難以使用該方法，而平面波成像為了提高橫向分辨率，通常需要采集多個(gè)偏轉(zhuǎn)角度的圖像進(jìn)行多角度空間復(fù)合。其次，在2個(gè)Kernel的方法中，Kernel1函數(shù)算出的延時(shí)表只能存放于全局存儲(chǔ)器中，而線程對(duì)全局存儲(chǔ)器的訪問(wèn)存在著幾百個(gè)時(shí)鐘周期的延遲，相比于共享存儲(chǔ)器、寄存器等，訪問(wèn)速度很慢，這在一定程度上也影響了并行波束合成的速度。

為了改善2個(gè)Kernel方法的不足，在改進(jìn)的方法中利用Kernel2里已經(jīng)分配的線程來(lái)直接計(jì)算延時(shí)，因?yàn)镚PU里分配的線程都是輕量級(jí)線程，因此單條線程的計(jì)算能力并不突出，如果將龐大的計(jì)算量分配給若干條線程來(lái)計(jì)算，會(huì)極大地拖慢程序的運(yùn)行速度。同時(shí)，需要為計(jì)算得到的延時(shí)值設(shè)計(jì)合理的存儲(chǔ)策略，讓其他的線程在進(jìn)行疊加計(jì)算時(shí)能夠達(dá)到較高的訪問(wèn)速度，保證1個(gè)Kernel的方法也能夠保持相同的計(jì)算能力。

利用1個(gè)Kernel實(shí)現(xiàn)延時(shí)-疊加波束合成，需要用到GPU中的共享存儲(chǔ)器，之所以在2個(gè)Kernel的方法中不使用共享存儲(chǔ)器，是因?yàn)镵ernel2在運(yùn)行過(guò)程中，一條線程訪問(wèn)延時(shí)表所在的存儲(chǔ)器并讀取相應(yīng)的延時(shí)值的次數(shù)只有一次，不存在一條線程對(duì)一個(gè)延時(shí)值重復(fù)讀取的情況，因此在2個(gè)Kernel的方法中使用共享存儲(chǔ)器并不會(huì)提高波束合成的速度。

使用1個(gè)Kernel實(shí)現(xiàn)延時(shí)-疊加計(jì)算的Grid和Block設(shè)計(jì)及并行任務(wù)劃分如圖5所示。

圖5 1個(gè)Kernel并行波束合成對(duì)應(yīng)的Grid和Block設(shè)計(jì)及并行任務(wù)劃分Fig.5 Diagrammatic sketch of grid and block designs and task assignments in beamforming with one kernel

在改進(jìn)的方法中舍棄了Kernel1，在Kernel2已經(jīng)分配的線程中分配一部分線程來(lái)計(jì)算延時(shí)值，并將計(jì)算得到的延時(shí)值存入共享存儲(chǔ)器中，避免了將延時(shí)值存放在訪問(wèn)速度較慢的全局存儲(chǔ)器中。在改進(jìn)Kernel2中，一個(gè)Block的任務(wù)是對(duì)同一個(gè)深度的N+1個(gè)數(shù)據(jù)進(jìn)行波束合成，并且同一個(gè)深度的N+1個(gè)數(shù)據(jù)在波束合成時(shí)共用2R+1個(gè)延時(shí)值。因此在劃分并行任務(wù)時(shí)，分配Kernel2的每一個(gè)Block里的2R+1條線程來(lái)并行地計(jì)算這個(gè)Block所對(duì)應(yīng)深度的2R+1個(gè)延時(shí)值，再并行地將這些延時(shí)值存放于所在共享存儲(chǔ)器中，經(jīng)過(guò)同步機(jī)制后，Block里的所有線程進(jìn)行疊加計(jì)算時(shí)，只需要從各自的共享存儲(chǔ)器中讀取相應(yīng)的延時(shí)值并使用即可。這樣的任務(wù)劃分相比于前者，線程訪問(wèn)讀取延時(shí)值的速度更快、效率更高。實(shí)驗(yàn)結(jié)果也證明，盡管Block里有一部分線程因?yàn)橛?jì)算延時(shí)值而增加了計(jì)算負(fù)擔(dān)，但因?yàn)槭褂昧嗽L問(wèn)速度更快(相比于全局存儲(chǔ)器)的共享存儲(chǔ)器，波束合成的速度并沒(méi)有降低。

1.3 實(shí)驗(yàn)設(shè)計(jì)

在仿體實(shí)驗(yàn)中，使用Verasonics V-1超聲數(shù)據(jù)采集系統(tǒng)和L11-4線陣探頭采集通道數(shù)據(jù)，探頭中心頻率7.5 MHz。實(shí)驗(yàn)中進(jìn)行對(duì)比的波束合成處理器分別是Intel公司的Xeon CPU(主頻3.10 GHz)和NVIDIA公司的Tesla K20c GPU(2 496個(gè)CUDA核，5 GB內(nèi)存)，采集通道RF數(shù)據(jù)的各項(xiàng)參數(shù)如表2所示。為了驗(yàn)證GPU波束合成結(jié)果的正確性，對(duì)合成的結(jié)果沒(méi)有進(jìn)行濾波、角度復(fù)合等優(yōu)化，因此獲得的圖像有較強(qiáng)的噪聲干擾，對(duì)比度和信噪比也不及一般的超聲圖像。其中，接收孔徑33代表合成每一個(gè)像素點(diǎn)時(shí)用到的通道數(shù)，接收孔徑大小為9.5 mm。

表2 線陣探頭采集通道RF參數(shù)及波束合成參數(shù)

1.4 評(píng)價(jià)方法

本研究中對(duì)于實(shí)驗(yàn)結(jié)果的評(píng)價(jià)方法如下：

1) 均方誤差。用歸一化的均方誤差來(lái)評(píng)估波束合成結(jié)果的正確性，其表達(dá)式如下：

(2)

式中：xi,j是參考數(shù)據(jù)集，在本研究中是采用傳統(tǒng)波束合成方法得到的圖像矩陣；yi,j是實(shí)驗(yàn)數(shù)據(jù)集，在本研究中是采用GPU波束合成方法得到的圖像矩陣；N是數(shù)據(jù)集的大??；mean(xi,j)表示求平均計(jì)算。

2) 計(jì)算時(shí)間。將超聲采集系統(tǒng)獲得的100幀通道數(shù)據(jù)進(jìn)行打包處理，再一次性送入波束合成器(CPU或GPU)中，計(jì)算波束合成器采用不同方法合成完畢100幀圖像的時(shí)間，以此作為衡量波束合成器計(jì)算性能的依據(jù)。

3) 算法復(fù)雜度和計(jì)算量。從算法本身的復(fù)雜程度對(duì)比，以及不同算法合成100幀圖像時(shí)所需要的計(jì)算量，對(duì)3種不同的方法進(jìn)行討論。

3 結(jié)果

仿體實(shí)驗(yàn)的結(jié)果如圖6所示，其中圖6(a)是采用普通串行算法計(jì)算得到的B模式圖像，圖6(b)是采用GPU并行算法中的2個(gè)Kernel方法計(jì)算得到的結(jié)果，圖6(c)是1個(gè)Kernel方法計(jì)算得到的結(jié)果。

1) 均方誤差。以圖6(a)作為參考數(shù)據(jù)集，將圖6(b)、(c)分別和參考數(shù)據(jù)集做均方誤差計(jì)算，得到的結(jié)果均為0.002 5%。

2) 計(jì)算時(shí)間。普通波束合成方法和兩種并行波束合成方法的計(jì)算能力對(duì)比如表3所示。可以看到，普通方法耗時(shí)4 500ms，幀頻為22f/s左右。2個(gè)Kernel的方法耗時(shí)為45.92ms，計(jì)算幀頻2 178f/s左右；使用1個(gè)Kernel的方法耗時(shí)40.76ms，幀頻達(dá)到2 453f/s左右。相比于普通方法，兩種并行合成方法分別加速99倍和111倍。相比2個(gè)Kernel的方法，1個(gè)Kernel的方法又進(jìn)一步提速12.68%。

表3 GPU并行波束合成和普通波束合成方法計(jì)算能力對(duì)比

4 討論

4.1 算法復(fù)雜度和計(jì)算量

基于CPU串行計(jì)算的波束合成方法在合成每一個(gè)像素點(diǎn)時(shí)，需要分別計(jì)算這些像素點(diǎn)的延時(shí)值，再進(jìn)行疊加計(jì)算。對(duì)于待合成圖像中的每一個(gè)像素點(diǎn)，需要循環(huán)重復(fù)上述過(guò)程。本研究提出的第2種GPU波束合成方法，即基于1個(gè)Kernel的方法，從算法復(fù)雜度的角度而言，和基于CPU串行計(jì)算的方法類似，即在合成每一個(gè)像素點(diǎn)時(shí)，該點(diǎn)的延時(shí)值是單獨(dú)分別計(jì)算和儲(chǔ)存的，不同的是后者將該方法進(jìn)行了并行化處理，并利用共享存儲(chǔ)器來(lái)提高存儲(chǔ)效率。綜合言之，這兩種方法的計(jì)算思路是一致的，但基于GPU的方法需要考慮并行分配線程和數(shù)據(jù)的存儲(chǔ)策略。

本研究中提出的基于2個(gè)Kernel的方法，需要先建立1個(gè)Kernel來(lái)計(jì)算延時(shí)表并儲(chǔ)存在全局存儲(chǔ)器中，再由第2個(gè)Kernel調(diào)用。從算法角度而言，該方法的復(fù)雜度相對(duì)較高，而且相比1個(gè)Kernel的方法，更多地利用了全局存儲(chǔ)器，因此在一定程度上降低了數(shù)據(jù)的讀取效率。

從數(shù)據(jù)量的角度考慮，基于CPU串行計(jì)算的方法和兩種基于GPU并行計(jì)算的方法，其數(shù)據(jù)的輸入量和輸出量是相同的，其中輸入是100幀的通道數(shù)據(jù)，輸出是100幀波束合成之后的射頻數(shù)據(jù)(數(shù)據(jù)格式也相同)。從中間生成的臨時(shí)數(shù)據(jù)而言，基于CPU的方法和基于GPU的1個(gè)Kernel的方法在合成每一個(gè)像素點(diǎn)時(shí)，只會(huì)產(chǎn)生對(duì)應(yīng)于該點(diǎn)的延時(shí)值數(shù)列，存儲(chǔ)于共享存儲(chǔ)器中并隨著波束合成的進(jìn)行而不斷更新。而基于GPU的2個(gè)Kernel的方法在計(jì)算過(guò)程中會(huì)產(chǎn)生一個(gè)延時(shí)表，它會(huì)根據(jù)圖像矩陣的大小而占據(jù)不同的存儲(chǔ)空間。

4.2 應(yīng)用靈活性

使用2個(gè)Kernel的方法需要先計(jì)算出延時(shí)表，是基于無(wú)角度偏轉(zhuǎn)發(fā)射/接收模式的平面波成像在同一深度處有完全相同的延時(shí)模式而采用的簡(jiǎn)便算法，它所適用的成像方法范圍較窄，因此當(dāng)遇到角度偏轉(zhuǎn)發(fā)射/接收模式的平面波成像或者是其他先進(jìn)成像方法時(shí)，這樣的方法使用并不方便。而1個(gè)Kernel的方法通過(guò)分配出部分線程直接計(jì)算所在的線程塊所需要的延時(shí)值，提高了應(yīng)用的靈活性。

4.3 創(chuàng)新和不足

上述兩種方法與已有的并行波束合成算法相比，不同之處在于，現(xiàn)有并行波束合成方法一般是同時(shí)對(duì)不同通道的RF數(shù)據(jù)進(jìn)行波束合成，而本研究所采用的方法考慮到了延時(shí)-疊加波束合成算法在同一深度上具有相同的延時(shí)模式，可以共享一組延時(shí)值，從而實(shí)現(xiàn)同一個(gè)深度上的并行波束合成。在大部分超聲圖像中，縱向上一個(gè)通道的采樣點(diǎn)數(shù)要遠(yuǎn)遠(yuǎn)大于橫向上的通道數(shù)，該方法并行地合成各個(gè)深度上的數(shù)據(jù)，相比現(xiàn)有的并行波束合成方法，具有更大的加速潛能。

在本研究中，根據(jù)延時(shí)值計(jì)算和存儲(chǔ)策略的不同，設(shè)計(jì)并實(shí)現(xiàn)了兩種基于GPU的并行波束合成方法，兩種方法最大的區(qū)別是對(duì)共享存儲(chǔ)器的使用。由于不同GPU的共享存儲(chǔ)器空間存在一定的差別(一般為48 KB)，因此采用第二種方法時(shí)需要考慮具體所使用的GPU共享存儲(chǔ)器大小，避免發(fā)生數(shù)據(jù)的溢出，導(dǎo)致計(jì)算結(jié)果錯(cuò)誤。

該項(xiàng)工作還有可以進(jìn)一步改進(jìn)和突破的地方，如進(jìn)一步優(yōu)化數(shù)據(jù)在GPU中的存儲(chǔ)策略，減少待波束合成數(shù)據(jù)因存放于全局存儲(chǔ)器中而造成的較高的訪問(wèn)延遲。另外，加入角度復(fù)合之后，數(shù)據(jù)量的成倍增加會(huì)導(dǎo)致計(jì)算幀頻的相應(yīng)降低，如何保證GPU波束合成器在數(shù)據(jù)量增加的情況下仍然保持較高的計(jì)算能力也是今后工作的重點(diǎn)之一。

5 結(jié)論

本研究基于延時(shí)-疊加波束合成算法的并行性分析，設(shè)計(jì)出了兩種基于GPU并行計(jì)算的超聲波束合成方法——基于2個(gè)Kernel和基于1個(gè)Kernel的并行波束合成方法，并分別編程實(shí)現(xiàn)。兩種方法均達(dá)到了較高的計(jì)算能力，相比普通的合成方法大幅度提高了計(jì)算幀頻。

[1] Bruneel C, Torguet R, Rouvaen KM, et al. Ultrafast echotomographic system using optical processing of ultrasonic signals[J]. Applied Physics Letters, 1977, 30(8): 371-373.

[2] Delannoy B, Torguet R, Bruneel C, et al. Acoustical image reconstruction in parallel‐processing analog electronic systems[J]. Journal of Applied Physics, 1979, 50(5): 3153-3159.

[3] Shattuck DP, Weinshenker MD, Smith SW, et al. Explososcan: A parallel processing technique for high speed ultrasound imaging with linear phased arrays[J]. The Journal of the Acoustical Society of America, 1984, 75(4): 1273-1282.

[4] Lu J, Greenleaf JF. Pulse-echo imaging using a nondiffracting beam transducer[J]. Ultrasound in Medicine & Biology, 1991, 17(3): 265-281.

[5] Tanter M, Fink M. Ultrafast imaging in biomedical ultrasound[J]. IEEE Trans Ultrasonics, Ferroelectrics, and Frequency Control, 2014, 61(1): 102-119.

[6] Souquet J, Bercoff J. Ultrafast ultrasound imaging[J]. Ultrasound in Medicine & Biology, 2011, 37(8): S17.

[7] Montaldo G, Tanter M, Bercoff J, et al. Coherent plane-wave compounding for very high frame rate ultrasonography and transient elastography[J]. IEEE Trans Ultrasonics, Ferroelectrics, and Frequency Control, 2009, 56(3): 489-506.

[8] 夏春蘭, 石丹, 劉東權(quán). 基于 CUDA 的超聲 B 模式成像[J]. 計(jì)算機(jī)應(yīng)用研究, 2011, 28(6): 2011-2015.

[9] So HKH, Chen J, Yiu BYS, et al. Medical ultrasound imaging: To GPU or not to GPU? [J]. IEEE Micro, 2011, 31(5): 54-65.

[10] Yiu BYS, Tsang IKH, Yu ACH. Real-time GPU-based software beamformer designed for advanced imaging methods research[C] // 2010 IEEE International Ultrasonics Symposium (IUS). San Diego: IEEE-UFFC, 2010: 1920-1923.

[11] Yiu BYS, Tsang IKH, Yu ACH. GPU-based beamformer: Fast realization of plane wave compounding and synthetic aperture imaging[J]. IEEE Trans Ultrasonics, Ferroelectrics, and Frequency Control, 2011, 58(8): 1698-1705.

[12] Chen J, Yiu BYS, So HKH, et al. Real-time GPU-based adaptive beamformer for high quality ultrasound imaging[C] // Clemens Ruppel, Ken-ya Hashimoto, eds. 2011 IEEE International Ultrasonics Symposium (IUS). Orlando: IEEE, 2011: 474-477.

[13] Cook S. CUDA Programming: A Developer′s Guide to Parallel Computing with GPUs [M]. Waltham: Elsevier, 2012.

[14] Farber R. CUDA Application Design and Development [M]. Waltham: Elsevier, 2011.

[15] Kirk DB, Wen-mei WH. Programming massively parallel processors: a hands-on approach [M]. Waltham: Elsevier, 2012.

[16] 張舒, 褚艷利. GPU 高性能運(yùn)算之 CUDA[M]. 北京: 中國(guó)水利水電出版社, 2009.

Ultrasound Beamforming Based on GPU Parallel Computation

Chen Yinran1,2He Qiong1,2Luo Jianwen1,2*

1(DepartmentofBiomedicalEngineering,SchoolofMedicine,TsinghuaUniversity,Beijing100084,China)2(CenterforBiomedicalImagingResearch,TsinghuaUniversity,Beijing100084,China)

The conventional line-by-line focused mode restricts the improvement of frame rate in ultrasound imaging. In plane wave imaging, each image is obtained by only one transmission and reception, and thus ultrafast imaging can be achieved. However, the current beamformers cannot meet the demand of ultrafast ultrasound imaging for massive computation. In this paper, the feasibility analysis of parallel computation in delay-and-sum beamforming algorithm was performed and two beamforming methods based on graphics processing unit (GPU) parallel computation including the 2 Kernel-based and 1 Kernel-based parallel beamformers were designed and implemented. The main differences between these methods were the calculations and data storage strategies of time delays in beamforming. Phantom experiments demonstrate that the computational frame rates of each method was 2,178 frames per second and 2,453 frames per second, respectively. Each of the methods obtained a speed up ratio of 99 and 111 compared to the normal method, which demonstrated that the GPU-based beamformer could significantly improve the calculating capability.

ultrasound imaging; beamforming; parallel computing; GPU

10.3969/j.issn.0258-8021. 2016. 06.006

2016-03-21， 錄用日期:2016-06-23

國(guó)家自然科學(xué)基金(61271131)；高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金(20130002110061)

R318

A

0258-8021(2016) 06-0677-07

*通信作者(Corresponding author)， E-mail: luo_jianwen@tsinghua.edu.cn