摘  要： 應用GPU通用高性能編程技術(shù)實現(xiàn)一種加速地震疊前時間偏移的新方法。該技術(shù)是地震勘探處理的常規(guī)流程，其核心算法具有計算密集、數(shù)據(jù)獨立性強、并行性高等特點。通過性能剖析獲得其計算熱點，通過CUDA技術(shù)對其進行并行化改造，并利用CUDA的流技術(shù)實現(xiàn)CPU到GPU的異步傳輸。通過集群環(huán)境下的性能測試，應用GPU并行化的PSTM程序可明顯縮短運行時間。
關(guān)鍵詞： 疊前時間偏移；地震勘探；圖像處理器；計算統(tǒng)一設備架構(gòu)；流；并行計算

1 PSTM技術(shù)
1.1 PSTM簡介
    PSTM是復雜構(gòu)造成像最有效的方法之一，能適應縱向速度變化較大的情況[1]，適用于大傾角的偏移成像[2]。PSTM經(jīng)過多年研究，20世紀90年代初期開始初步應用，中后期在不少探區(qū)的地震勘探中發(fā)揮了重要作用，進入21世紀后開始了較為廣泛的應用，目前部分處理公司和計算中心已把該技術(shù)作為常規(guī)軟件加入到常規(guī)處理流程中，成為獲取保幅信息實現(xiàn)屬性分析、AVO/AVA/AVP反演和其他參數(shù)反演的重要步驟和依據(jù)。
    PSTM方法主要分為兩類，即用于準確構(gòu)造成像的PSTM和振幅保持PSTM，每一類方法都有兩種實現(xiàn)方式：Kirchhoff型和波動方程型，而目前工業(yè)界最成熟的是Kirchhoff型PSTM[3]。PSTM每輸出一個地震道，就是一次海量運算。以1 ms采樣、6 s數(shù)據(jù)為例，一個地震道的輸出需要至少1 000萬道甚至更多（偏移孔徑?jīng)Q定）的輸入道，每一個點要做兩次均方根運算以及兩次加法運算，振幅補償兩次乘法運算。如此計算下來，實現(xiàn)一道偏移需要1 000 000×6 000×2×(平方+加法+乘法)次數(shù)學運算，計算量和需要處理的數(shù)據(jù)量都極其巨大[4，5]。目前，人們往往使用大規(guī)模的服務器集群來進行疊前偏移處理，其原理是將數(shù)據(jù)先分配到各個CPU核上，然后由各個CPU核單獨進行計算，最后將結(jié)果匯總輸出。這種做法消耗了大量的時間、電力和維護費用。而且，隨著人們對石油勘探地震資料處理的周期要求越來越短，精度要求越來越高，服務器集群的規(guī)模越做越大，在系統(tǒng)構(gòu)建成本、數(shù)據(jù)中心機房空間、內(nèi)存和I/O帶寬、功耗散熱和電力限制、可管理性、編程簡易性、擴展性、管理維護費用等方面都面臨著巨大的挑戰(zhàn)。
1.2 PSTM研究現(xiàn)狀
    自從20世紀90年代以來，疊前時間偏移在國外取得了很大發(fā)展。在理論研究方面，Bleistein、Bortfeld和Hubral等進行了一系列有關(guān)真振幅疊前時間偏移理論的研究工作[6]，Schneider給出了Kirchhoff型真振幅偏移權(quán)函數(shù)的一般公式。Graham A.Winbow（1999）推出了控制振幅的三維疊前時間偏移的權(quán)函數(shù)的顯式公式，并且利用真振幅權(quán)函數(shù)估計進行了振幅補償。另外，在權(quán)函數(shù)改進的基礎上，提出了高精度的彎曲射線法繞射走時計算方法。在應用方面，最近幾年，在常規(guī)疊前時間偏移基礎上，研究開發(fā)了多種保幅型疊前時間偏移軟件，尤其是Kirchhoff保幅型疊前時間偏移軟件取得了巨大成功。隨著GPU的出現(xiàn)，利用GPU的多核計算資源，實現(xiàn)對PSTM的并行處理，加速PSTM的執(zhí)行時間，成為一個研究熱門，國內(nèi)外多家研究及工程機構(gòu)皆針對此技術(shù)進行了研究。
1.3 GPU通用高性能編程技術(shù)
    多核技術(shù)已經(jīng)代替單純CPU頻率的提升而成為計算機性能發(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更高的晶體管密度和計算能力，相較FPGA等的實現(xiàn)，GPU是一種非常成熟的商業(yè)芯片，有利于應用的快速部署，并表現(xiàn)出更好的性價比。從2000年左右開始，有部分研究者開始采用GPU強大的圖像處理能力來進行通用計算[7]。傳統(tǒng)的GPGPU(General-Purpose Computation on Graphics Processing Units)的處理方法是把通用的計算問題轉(zhuǎn)換為GPU中能處理的圖形問題，然后通過圖形編程語言來進行編程。這樣的編程模式使得開發(fā)難度很高，開發(fā)者必須掌握強大的圖形處理能力。
    從2006年開始，NVIDIA開始研究通用的GPU處理架構(gòu)并推出計算統(tǒng)一設備架構(gòu)CUDA(Compute Unified Device Archetecture)，使得開發(fā)者不需要有很強的圖形開發(fā)能力，而是采用擴展了的C語言對GPU進行直接編程。NVIDIA把頂點處理和面處理整合到同一個處理芯片上，從而獲得GPU硬件在架構(gòu)上的通用性。經(jīng)過幾年的發(fā)展，CUDA架構(gòu)更加靈活，更適合通用計算。這使得CUDA方案成為異構(gòu)并行模型中的佼佼者。該方案入門門檻低，程序移植效率高，加速比好。圖1描述了CUDA技術(shù)體系結(jié)構(gòu)。

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

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

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

    既然PSTM_Kernel為PSTM的熱點計算部分，如果對PSTM_Kernel進行并行化改造，它的運行時間將大大減少，則整個PSTM的運行時間將也隨之明顯減少，其性能將大大提高，所以在本文的研究中，將對PSTM_Kernel進行并行化的改造，以加速PSTM的運行效率。
3 PSTM_Kernel并行化改造
3.1 PSTM_Kernel串行算法分析
    PSTM_Kernel串行算法的基本思想是：對于每一道地震道數(shù)據(jù)而言，PSTM_Kernel都循環(huán)計算成像空間中每一個點的走時，即每一個成像點到此道地震道數(shù)據(jù)對應的炮點的走時，和每一個成像點到此道地震道數(shù)據(jù)對應的接收點的走時[7]。然后根據(jù)這兩個走時之和取出對應的地震道數(shù)據(jù)，更新此成像點的數(shù)據(jù)，實現(xiàn)疊前時間偏移計算。
    PSTM_Kernel串行算法的核心是按照成像空間X軸方向、Y軸方向、Z軸方向三層循環(huán)進行實現(xiàn)的，其具體過程如下：
    (1)循環(huán)取出每個成像點在成像空間中對應的X、Y、Z坐標，確定成像點位置；
    (2)根據(jù)成像點位置計算出其到炮點和接收點的走時；
    (3)根據(jù)這兩個走時之和取出對應的地震道數(shù)據(jù)，然后更新此成像點在成像空間中的數(shù)據(jù)，實現(xiàn)疊前時間偏移計算。
    從以上串行算法分析，每一個成像點到炮點和接收點的走時計算是數(shù)據(jù)獨立的，并不存在依賴性，因此存在較好的并行性，適合GPU等細粒度并行體系結(jié)構(gòu)的加速。在此采用CUDA技術(shù)進行實現(xiàn)。
3.2 CUDA并行實現(xiàn)
3.2.1 線程計算粒度選擇
    應用CUDA模型加速通用計算，涉及線程模型及對應的計算粒度問題。如果線程計算的粒度太小，例如成像空間的每一個點對應一個線程進行計算，則線程總數(shù)為成像空間四個維度數(shù)值的乘積，即NX×NY×NOFF×NZ。由于每道道數(shù)據(jù)并非對成像空間的所有點都有貢獻，因此細粒度的并行會產(chǎn)生較多無效計算線程。這種情況下，線程計算粒度太小反而會導致線程的并發(fā)度不高。如果線程計算粒度太大，一個線程進行多點的循環(huán)計算，則線程數(shù)過少，不能使流處理器滿負荷運行，性能同樣不能達到最優(yōu)。通過性能測試發(fā)現(xiàn)，最佳計算粒度為：32個線程負責計算NZ個點，即32個線程負責計算成像空間Z方向上一條線上的所有點，此時計算性能最為理想。
3.2.2 線程模型選擇
    線程模型是用來明確CUDA程序內(nèi)核的執(zhí)行配置。根據(jù)GPU硬件資源，定義網(wǎng)格和線程塊，好的線程模型能大大提升程序的性能。
    (1)網(wǎng)格（grid）的定義：即把所有線程如何進行分塊，定義線程塊數(shù)和線程塊的組織方式，這里根據(jù)成像空間的X-Y平面來劃分線程塊，其具體定義為dimGrid(NX，NY/10)；
    (2)線程塊（block）的定義：即定義一個線程塊有多少個線程和線程的組織方式，根據(jù)內(nèi)核所需的寄存器數(shù)和共享內(nèi)存數(shù)量，定義一個線程塊為320個線程，其具體定義為dimBlock(32，10)；
    (3)線程模型描述：grid和block都是定義為二維的，線程總數(shù)為NX×NY×32，線程塊數(shù)為(NX×NY)/10。每個block的320個線程處理10個(x，y)坐標所對應的Z軸的所有點，即處理Z軸的10條線。如果把每個block的320個線程按照32個線程進行分組，每一組為一個warp，則整個block有10個warp，第0個block的第0個warp，即threadIdx.y等于0的線程處理第0個(x，y)坐標對應的Z軸的第一條線；第0個block的第1個warp，即threadIdx.y等于1的線程處理第1個(x，y)坐標對應的Z軸的另一條線，依此類推，第0個block的第9個warp，處理第9個(x，y)坐標對應的Z軸的一條線。同理，其他block處理另外10個(x，y)坐標對應的Z軸的10條線。10條線并行進行走時計算，可以使計算更均衡，性能更高。
3.3 CPU與GPU的異步IO實現(xiàn)
    基于上述的CUDA并行化改造，顯著提升了計算性能。而隨之帶來的問題是GPU加速所必然造成的IO傳輸問題。PSTM_Kernel執(zhí)行之前，需要FFT計算模塊對輸入道進行預處理，之后這些數(shù)據(jù)需要從系統(tǒng)內(nèi)存?zhèn)鬟f至GPU顯存。這是GPU計算的引入所帶來的額外IO開銷。如果采用同步方式進行數(shù)據(jù)傳輸，即：對一道地震道數(shù)據(jù)而言，完成一道FFT的計算，就將數(shù)據(jù)傳輸給GPU再進行疊前時間偏移計算，那么IO開銷將會很大，對性能造成很大影響。
    本文采用CPU與GPU的異步傳輸方式來試圖減少上述PSTM執(zhí)行的總時間。采用CUDA的流(stream)技術(shù)，利用雙流雙緩沖，實現(xiàn)CPU與GPU的異步IO傳輸，其邏輯結(jié)構(gòu)圖如圖２所示。

4 地震資料測試
4.1 測試環(huán)境及數(shù)據(jù)
    通過一套小規(guī)模集群系統(tǒng)測試GPU加速PSTM的效果。具體各項參數(shù)如表3所示。

4.2 性能結(jié)果
    為了保證測試性能結(jié)果的穩(wěn)定性，對上述體偏作業(yè)進行了3次測試，CPU多線程版PSTM在集群上運行3次的平均時間是110 908 s，而GPU版PSTM在集群上運行上述同樣體偏作業(yè)3次的平均時間是21 454 s，GPU版PSTM運行的性能是CPU多線程版PSTM的110 908/21 454=5倍。
4.3 成像效果圖比對
    通過上述作業(yè)，從獲得的效果圖來看，圖像不存在明顯差異，證明GPU加速PSTM運行結(jié)果的正確性。
    通過對PSTM算法的性能剖析，分析算法的熱點函數(shù)，并對該熱點函數(shù)進行GPU并行化改造。初步改造移植結(jié)果說明，使用GPU進行加速可顯著提高疊前時間偏移計算速度，實驗數(shù)據(jù)證明了加速效果。測試數(shù)據(jù)表明，對于一個完整的PSTM體偏作業(yè)，一個GPU節(jié)點的運行時間是一個采用最新雙路CPU節(jié)點，并運行商用多線程CPU版PSTM時間的1/5。由此可見，PSTM高并行度、無數(shù)據(jù)依賴性等特征使得其較適合GPU類設備的加速。
參考文獻
[1] 張占杰.疊前時間偏移技術(shù)及其在東海復雜地震資料處理中的應用[J].海洋石油，2007，27(3)：22-26.
[2] 王余慶.疊前偏移技術(shù)探討及應用[J].西北油氣勘探，2006，18(2)：31-39.
[3] 王棣.疊前時間偏移方法綜述[J].勘探地球物理進展，2004，27(5)：313-320.
[4] 羅銀河.Kirchhoff彎曲射線疊前時間偏移及應用[J].天然氣工業(yè)，2005，25(8)：35-37.
[5] 洪釗峰.GPU計算：在石油勘探領(lǐng)域的革命[EB/OL].[2009-5-13].http://server.it168.com/a2009/0513/276/
000000276186.shtml.
[6] 張麗艷.相對振幅保持的轉(zhuǎn)換波疊前時間偏移方法研究[J].石油地球物理勘探，2008，43(2)：30-36.
[7] 李博.地震疊前時間偏移的一種圖形處理器提速實現(xiàn)方法[J].地球物理學報，2009，52(1)：245-252.