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

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

(1.浪潮集團 高效能服務器和存儲技術國家重點實驗室，北京100085；2.香港浸會大學 計算機系，香港；3.中國石油集團 東方地球物理勘探有限責任公司，河北 涿州072751)

1 PSTM技術

1.1 PSTM簡介

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

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

1.2 PSTM研究現狀

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

1.3 GPU通用高性能編程技術

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

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

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

2 PSTM熱點分析

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

2.1 測試環境及數據

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

2.2 分析結果

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

表1 測試環境和測試數據

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

3 PSTM_Kernel并行化改造

3.1 PSTM_Kernel串行算法分析

PSTM_Kernel串行算法的基本思想是：對于每一道地震道數據而言，PSTM_Kernel都循環計算成像空間中每一個點的走時，即每一個成像點到此道地震道數據對應的炮點的走時，和每一個成像點到此道地震道數據對應的接收點的走時[7]。然后根據這兩個走時之和取出對應的地震道數據，更新此成像點的數據，實現疊前時間偏移計算。

PSTM_Kernel串行算法的核心是按照成像空間X軸方向、Y軸方向、Z軸方向三層循環進行實現的，其具體過程如下：

(1)循環取出每個成像點在成像空間中對應的 X、Y、Z坐標，確定成像點位置；

(2)根據成像點位置計算出其到炮點和接收點的走時；

(3)根據這兩個走時之和取出對應的地震道數據，然后更新此成像點在成像空間中的數據，實現疊前時間偏移計算。

從以上串行算法分析，每一個成像點到炮點和接收點的走時計算是數據獨立的，并不存在依賴性，因此存在較好的并行性，適合GPU等細粒度并行體系結構的加速。在此采用CUDA技術進行實現。

3.2 CUDA并行實現

3.2.1線程計算粒度選擇

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

3.2.2線程模型選擇

線程模型是用來明確CUDA程序內核的執行配置。根據GPU硬件資源，定義網格和線程塊，好的線程模型能大大提升程序的性能。

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

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

(3)線程模型描述：grid和block都是定義為二維的，線程總數為 NX×NY×32，線程塊數為(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條線并行進行走時計算，可以使計算更均衡，性能更高。

表2 PSTM軟件各部分運行時間表

圖2 異步IO邏輯結構圖

3.3 CPU與GPU的異步IO實現

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

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

4 地震資料測試

4.1 測試環境及數據

通過一套小規模集群系統測試GPU加速PSTM的效果。具體各項參數如表3所示。

表3 測試環境及測試數據

4.2 性能結果

為了保證測試性能結果的穩定性，對上述體偏作業進行了3次測試，CPU多線程版PSTM在集群上運行3次的平均時間是110 908 s，而GPU版PSTM在集群上運行上述同樣體偏作業3次的平均時間是21 454 s，GPU版PSTM運行的性能是CPU多線程版PSTM的110 908/21 454=5倍。

4.3 成像效果圖比對

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

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

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

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

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

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

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

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

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