一種使用GPU加速地震疊前時(shí)間偏移的方法

張 清 ，謝海波 ，趙開勇 ，，吳 慶 ，陳 維 ，王獅虎 ，遲旭光 ，褚曉文

(1.浪潮集團(tuán) 高效能服務(wù)器和存儲技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京100085；2.香港浸會大學(xué) 計(jì)算機(jī)系，香港；3.中國石油集團(tuán) 東方地球物理勘探有限責(zé)任公司，河北 涿州072751)

1 PSTM技術(shù)

1.1 PSTM簡介

PSTM是復(fù)雜構(gòu)造成像最有效的方法之一，能適應(yīng)縱向速度變化較大的情況[1]，適用于大傾角的偏移成像[2]。PSTM經(jīng)過多年研究，20世紀(jì) 90年代初期開始初步應(yīng)用，中后期在不少探區(qū)的地震勘探中發(fā)揮了重要作用，進(jìn)入21世紀(jì)后開始了較為廣泛的應(yīng)用，目前部分處理公司和計(jì)算中心已把該技術(shù)作為常規(guī)軟件加入到常規(guī)處理流程中，成為獲取保幅信息實(shí)現(xiàn)屬性分析、AVO/AVA/AVP反演和其他參數(shù)反演的重要步驟和依據(jù)。

PSTM方法主要分為兩類，即用于準(zhǔn)確構(gòu)造成像的PSTM和振幅保持PSTM，每一類方法都有兩種實(shí)現(xiàn)方式：Kirchhoff型和波動方程型，而目前工業(yè)界最成熟的是 Kirchhoff型 PSTM[3]。PSTM每輸出一個(gè)地震道，就是一次海量運(yùn)算。以1 ms采樣、6 s數(shù)據(jù)為例，一個(gè)地震道的輸出需要至少1 000萬道甚至更多（偏移孔徑?jīng)Q定）的輸入道，每一個(gè)點(diǎn)要做兩次均方根運(yùn)算以及兩次加法運(yùn)算，振幅補(bǔ)償兩次乘法運(yùn)算。如此計(jì)算下來，實(shí)現(xiàn)一道偏移需要 1 000 000×6 000×2×(平方+加法+乘法)次數(shù)學(xué)運(yùn)算，計(jì)算量和需要處理的數(shù)據(jù)量都極其巨大[4，5]。 目前，人們往往使用大規(guī)模的服務(wù)器集群來進(jìn)行疊前偏移處理，其原理是將數(shù)據(jù)先分配到各個(gè)CPU核上，然后由各個(gè)CPU核單獨(dú)進(jìn)行計(jì)算，最后將結(jié)果匯總輸出。這種做法消耗了大量的時(shí)間、電力和維護(hù)費(fèi)用。而且，隨著人們對石油勘探地震資料處理的周期要求越來越短，精度要求越來越高，服務(wù)器集群的規(guī)模越做越大，在系統(tǒng)構(gòu)建成本、數(shù)據(jù)中心機(jī)房空間、內(nèi)存和I/O帶寬、功耗散熱和電力限制、可管理性、編程簡易性、擴(kuò)展性、管理維護(hù)費(fèi)用等方面都面臨著巨大的挑戰(zhàn)。

1.2 PSTM研究現(xiàn)狀

自從20世紀(jì)90年代以來，疊前時(shí)間偏移在國外取得了很大發(fā)展。在理論研究方面，Bleistein、Bortfeld和Hubral等進(jìn)行了一系列有關(guān)真振幅疊前時(shí)間偏移理論的研究工作[6]，Schneider給出了 Kirchhoff型真振幅偏移權(quán)函數(shù)的一般公式。Graham A.Winbow（1999）推出了控制振幅的三維疊前時(shí)間偏移的權(quán)函數(shù)的顯式公式，并且利用真振幅權(quán)函數(shù)估計(jì)進(jìn)行了振幅補(bǔ)償。另外，在權(quán)函數(shù)改進(jìn)的基礎(chǔ)上，提出了高精度的彎曲射線法繞射走時(shí)計(jì)算方法。在應(yīng)用方面，最近幾年，在常規(guī)疊前時(shí)間偏移基礎(chǔ)上，研究開發(fā)了多種保幅型疊前時(shí)間偏移軟件，尤其是Kirchhoff保幅型疊前時(shí)間偏移軟件取得了巨大成功。隨著GPU的出現(xiàn)，利用GPU的多核計(jì)算資源，實(shí)現(xiàn)對PSTM的并行處理，加速PSTM的執(zhí)行時(shí)間，成為一個(gè)研究熱門，國內(nèi)外多家研究及工程機(jī)構(gòu)皆針對此技術(shù)進(jìn)行了研究。

1.3 GPU通用高性能編程技術(shù)

多核技術(shù)已經(jīng)代替單純CPU頻率的提升而成為計(jì)算機(jī)性能發(fā)展的主流。除了CPU技術(shù)之外，采用GPU、FPGA加速并與CPU異構(gòu)并行的方式在這兩年大放異彩。使用這種架構(gòu)的技術(shù)如ClearSpeed的FPGA加速方案、AMD的Stream方案、Mercury的Cell BE加速方案和NVIDIA的CUDA方案。

GPU已經(jīng)體現(xiàn)了較CPU更高的晶體管密度和計(jì)算能力，相較FPGA等的實(shí)現(xiàn)，GPU是一種非常成熟的商業(yè)芯片，有利于應(yīng)用的快速部署，并表現(xiàn)出更好的性價(jià)比。從2000年左右開始，有部分研究者開始采用GPU強(qiáng)大的圖像處理能力來進(jìn)行通用計(jì)算[7]。傳統(tǒng)的GPGPU(General-Purpose Computation on Graphics Processing Units)的處理方法是把通用的計(jì)算問題轉(zhuǎn)換為GPU中能處理的圖形問題，然后通過圖形編程語言來進(jìn)行編程。這樣的編程模式使得開發(fā)難度很高，開發(fā)者必須掌握強(qiáng)大的圖形處理能力。

從2006年開始，NVIDIA開始研究通用的GPU處理架構(gòu)并推出計(jì)算統(tǒng)一設(shè)備架構(gòu)CUDA(Compute Unified Device Archetecture)，使得開發(fā)者不需要有很強(qiáng)的圖形開發(fā)能力，而是采用擴(kuò)展了的C語言對GPU進(jìn)行直接編程。NVIDIA把頂點(diǎn)處理和面處理整合到同一個(gè)處理芯片上，從而獲得GPU硬件在架構(gòu)上的通用性。經(jīng)過幾年的發(fā)展，CUDA架構(gòu)更加靈活，更適合通用計(jì)算。這使得CUDA方案成為異構(gòu)并行模型中的佼佼者。該方案入門門檻低，程序移植效率高，加速比好。圖1描述了CUDA技術(shù)體系結(jié)構(gòu)。

2 PSTM熱點(diǎn)分析

通過分析PSTM，不難發(fā)現(xiàn)整個(gè)程序分為兩個(gè)部分，一是FFT計(jì)算部分（以下簡稱 FFT），二是 PSTM計(jì)算部分(以下簡稱 PSTM_Kernel)。FFT為 PSTM_Kernel作數(shù)據(jù)準(zhǔn)備，它主要是對輸入地震道數(shù)據(jù)進(jìn)行FFT計(jì)算，在這里只對PSTM的實(shí)際運(yùn)行情況進(jìn)行定量分析，找到PSTM中的熱點(diǎn)計(jì)算部分。

2.1 測試環(huán)境及數(shù)據(jù)

測試環(huán)境包括硬件環(huán)境、軟件環(huán)境、運(yùn)行軟件；測試數(shù)據(jù)為輸入地震道數(shù)據(jù)集；對于成像空間而言，它是四維的，其中第一維為X軸方向的大小NX，第二維為Y軸方向的大小NY，第三維為炮點(diǎn)與接收點(diǎn)的偏移范圍大小NOFF，第四維為Z軸方向的大小NZ。為準(zhǔn)確測試出熱點(diǎn)，特選擇一個(gè)線偏成像空間和一個(gè)體偏成像空間進(jìn)行熱點(diǎn)測試，具體各項(xiàng)參數(shù)如表1所示。

2.2 分析結(jié)果

為了保證數(shù)據(jù)的準(zhǔn)確性，對線偏和體偏分別進(jìn)行10次測試，計(jì)算10次PSTM運(yùn)行的平均CPU使用時(shí)間和它的兩個(gè)計(jì)算部分平均執(zhí)行時(shí)間，測試結(jié)果如表2所示。不難發(fā)現(xiàn)，PSTM_Kernel是整個(gè)PSTM的核心計(jì)算部分，占到整個(gè)運(yùn)行時(shí)間的90%以上，所以PSTM_Kernel為PSTM的熱點(diǎn)計(jì)算部分。

表1 測試環(huán)境和測試數(shù)據(jù)

既然PSTM_Kernel為PSTM的熱點(diǎn)計(jì)算部分，如果對PSTM_Kernel進(jìn)行并行化改造，它的運(yùn)行時(shí)間將大大減少，則整個(gè)PSTM的運(yùn)行時(shí)間將也隨之明顯減少，其性能將大大提高，所以在本文的研究中，將對PSTM_Kernel進(jìn)行并行化的改造，以加速PSTM的運(yùn)行效率。

3 PSTM_Kernel并行化改造

3.1 PSTM_Kernel串行算法分析

PSTM_Kernel串行算法的基本思想是：對于每一道地震道數(shù)據(jù)而言，PSTM_Kernel都循環(huán)計(jì)算成像空間中每一個(gè)點(diǎn)的走時(shí)，即每一個(gè)成像點(diǎn)到此道地震道數(shù)據(jù)對應(yīng)的炮點(diǎn)的走時(shí)，和每一個(gè)成像點(diǎn)到此道地震道數(shù)據(jù)對應(yīng)的接收點(diǎn)的走時(shí)[7]。然后根據(jù)這兩個(gè)走時(shí)之和取出對應(yīng)的地震道數(shù)據(jù)，更新此成像點(diǎn)的數(shù)據(jù)，實(shí)現(xiàn)疊前時(shí)間偏移計(jì)算。

PSTM_Kernel串行算法的核心是按照成像空間X軸方向、Y軸方向、Z軸方向三層循環(huán)進(jìn)行實(shí)現(xiàn)的，其具體過程如下：

(1)循環(huán)取出每個(gè)成像點(diǎn)在成像空間中對應(yīng)的 X、Y、Z坐標(biāo)，確定成像點(diǎn)位置；

(2)根據(jù)成像點(diǎn)位置計(jì)算出其到炮點(diǎn)和接收點(diǎn)的走時(shí)；

(3)根據(jù)這兩個(gè)走時(shí)之和取出對應(yīng)的地震道數(shù)據(jù)，然后更新此成像點(diǎn)在成像空間中的數(shù)據(jù)，實(shí)現(xiàn)疊前時(shí)間偏移計(jì)算。

從以上串行算法分析，每一個(gè)成像點(diǎn)到炮點(diǎn)和接收點(diǎn)的走時(shí)計(jì)算是數(shù)據(jù)獨(dú)立的，并不存在依賴性，因此存在較好的并行性，適合GPU等細(xì)粒度并行體系結(jié)構(gòu)的加速。在此采用CUDA技術(shù)進(jìn)行實(shí)現(xiàn)。

3.2 CUDA并行實(shí)現(xiàn)

3.2.1線程計(jì)算粒度選擇

應(yīng)用CUDA模型加速通用計(jì)算，涉及線程模型及對應(yīng)的計(jì)算粒度問題。如果線程計(jì)算的粒度太小，例如成像空間的每一個(gè)點(diǎn)對應(yīng)一個(gè)線程進(jìn)行計(jì)算，則線程總數(shù)為成像空間四個(gè)維度數(shù)值的乘積，即NX×NY×NOFF×NZ。由于每道道數(shù)據(jù)并非對成像空間的所有點(diǎn)都有貢獻(xiàn)，因此細(xì)粒度的并行會產(chǎn)生較多無效計(jì)算線程。這種情況下，線程計(jì)算粒度太小反而會導(dǎo)致線程的并發(fā)度不高。如果線程計(jì)算粒度太大，一個(gè)線程進(jìn)行多點(diǎn)的循環(huán)計(jì)算，則線程數(shù)過少，不能使流處理器滿負(fù)荷運(yùn)行，性能同樣不能達(dá)到最優(yōu)。通過性能測試發(fā)現(xiàn)，最佳計(jì)算粒度為：32個(gè)線程負(fù)責(zé)計(jì)算NZ個(gè)點(diǎn)，即32個(gè)線程負(fù)責(zé)計(jì)算成像空間Z方向上一條線上的所有點(diǎn)，此時(shí)計(jì)算性能最為理想。

3.2.2線程模型選擇

線程模型是用來明確CUDA程序內(nèi)核的執(zhí)行配置。根據(jù)GPU硬件資源，定義網(wǎng)格和線程塊，好的線程模型能大大提升程序的性能。

(1)網(wǎng)格（grid）的定義：即把所有線程如何進(jìn)行分塊，定義線程塊數(shù)和線程塊的組織方式，這里根據(jù)成像空間的X-Y平面來劃分線程塊，其具體定義為dimGrid(NX，NY/10)；

(2)線程塊（block）的定義：即定義一個(gè)線程塊有多少個(gè)線程和線程的組織方式，根據(jù)內(nèi)核所需的寄存器數(shù)和共享內(nèi)存數(shù)量，定義一個(gè)線程塊為320個(gè)線程，其具體定義為 dimBlock(32，10)；

(3)線程模型描述：grid和block都是定義為二維的，線程總數(shù)為 NX×NY×32，線程塊數(shù)為(NX×NY)/10。 每個(gè)block的 320個(gè)線程處理 10個(gè)(x，y)坐標(biāo)所對應(yīng)的 Z軸的所有點(diǎn)，即處理Z軸的10條線。如果把每個(gè)block的320個(gè)線程按照32個(gè)線程進(jìn)行分組，每一組為一個(gè)warp，則整個(gè) block有 10個(gè) warp，第 0個(gè) block的第 0個(gè)warp，即 threadIdx.y等于 0的線程處理第 0個(gè)(x，y)坐標(biāo)對應(yīng)的Z軸的第一條線；第0個(gè)block的第1個(gè)warp，即threadIdx.y等于 1的線程處理第 1個(gè)(x，y)坐標(biāo)對應(yīng)的 Z軸的另一條線，依此類推，第0個(gè)block的第9個(gè)warp，處理第9個(gè)(x，y)坐標(biāo)對應(yīng)的Z軸的一條線。同理，其他block處理另外10個(gè)(x，y)坐標(biāo)對應(yīng)的Z軸的10條線。10條線并行進(jìn)行走時(shí)計(jì)算，可以使計(jì)算更均衡，性能更高。

表2 PSTM軟件各部分運(yùn)行時(shí)間表

圖2 異步IO邏輯結(jié)構(gòu)圖

3.3 CPU與GPU的異步IO實(shí)現(xiàn)

基于上述的CUDA并行化改造，顯著提升了計(jì)算性能。而隨之帶來的問題是GPU加速所必然造成的IO傳輸問題。PSTM_Kernel執(zhí)行之前，需要FFT計(jì)算模塊對輸入道進(jìn)行預(yù)處理，之后這些數(shù)據(jù)需要從系統(tǒng)內(nèi)存?zhèn)鬟f至GPU顯存。這是GPU計(jì)算的引入所帶來的額外IO開銷。如果采用同步方式進(jìn)行數(shù)據(jù)傳輸，即：對一道地震道數(shù)據(jù)而言，完成一道FFT的計(jì)算，就將數(shù)據(jù)傳輸給GPU再進(jìn)行疊前時(shí)間偏移計(jì)算，那么IO開銷將會很大，對性能造成很大影響。

本文采用CPU與GPU的異步傳輸方式來試圖減少上述PSTM執(zhí)行的總時(shí)間。采用CUDA的流(stream)技術(shù)，利用雙流雙緩沖，實(shí)現(xiàn)CPU與GPU的異步IO傳輸，其邏輯結(jié)構(gòu)圖如圖2所示。

4 地震資料測試

4.1 測試環(huán)境及數(shù)據(jù)

通過一套小規(guī)模集群系統(tǒng)測試GPU加速PSTM的效果。具體各項(xiàng)參數(shù)如表3所示。

表3 測試環(huán)境及測試數(shù)據(jù)

4.2 性能結(jié)果

為了保證測試性能結(jié)果的穩(wěn)定性，對上述體偏作業(yè)進(jìn)行了3次測試，CPU多線程版PSTM在集群上運(yùn)行3次的平均時(shí)間是110 908 s，而GPU版PSTM在集群上運(yùn)行上述同樣體偏作業(yè)3次的平均時(shí)間是21 454 s，GPU版PSTM運(yùn)行的性能是CPU多線程版PSTM的110 908/21 454=5倍。

4.3 成像效果圖比對

通過上述作業(yè)，從獲得的效果圖來看，圖像不存在明顯差異，證明GPU加速PSTM運(yùn)行結(jié)果的正確性。

通過對PSTM算法的性能剖析，分析算法的熱點(diǎn)函數(shù)，并對該熱點(diǎn)函數(shù)進(jìn)行GPU并行化改造。初步改造移植結(jié)果說明，使用GPU進(jìn)行加速可顯著提高疊前時(shí)間偏移計(jì)算速度，實(shí)驗(yàn)數(shù)據(jù)證明了加速效果。測試數(shù)據(jù)表明，對于一個(gè)完整的PSTM體偏作業(yè)，一個(gè)GPU節(jié)點(diǎn)的運(yùn)行時(shí)間是一個(gè)采用最新雙路CPU節(jié)點(diǎn)，并運(yùn)行商用多線程CPU版PSTM時(shí)間的1/5。由此可見，PSTM高并行度、無數(shù)據(jù)依賴性等特征使得其較適合GPU類設(shè)備的加速。

[1]張占杰.疊前時(shí)間偏移技術(shù)及其在東海復(fù)雜地震資料處理中的應(yīng)用[J].海洋石油，2007，27(3)：22-26.

[2]王余慶.疊前偏移技術(shù)探討及應(yīng)用[J].西北油氣勘探，2006，18(2)：31-39.

[3]王棣.疊前時(shí)間偏移方法綜述[J].勘探地球物理進(jìn)展，2004，27(5)：313-320.

[4]羅銀河.Kirchhoff彎曲射線疊前時(shí)間偏移及應(yīng)用[J].天然氣工業(yè)，2005，25(8)：35-37.

[5]洪釗峰.GPU計(jì)算：在石油勘探領(lǐng)域的革命[EB/OL].[2009-5-13].http://server.it168.com/a2009/0513/276/000000276186.shtml.

[6]張麗艷.相對振幅保持的轉(zhuǎn)換波疊前時(shí)間偏移方法研究[J].石油地球物理勘探，2008，43(2)：30-36.

[7]李博.地震疊前時(shí)間偏移的一種圖形處理器提速實(shí)現(xiàn)方法[J].地球物理學(xué)報(bào)，2009，52(1)：245-252.