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

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

1 PSTM技術

1.1 PSTM簡介

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

PSTM方法主要分為兩類，即用于準確構造成像的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核單獨進行計算，最后將結果匯總輸出。這種做法消耗了大量的時間、電力和維護費用。而且，隨著人們對石油勘探地震資料處理的周期要求越來越短，精度要求越來越高，服務器集群的規(guī)模越做越大，在系統(tǒng)構建成本、數(shù)據(jù)中心機房空間、內存和I/O帶寬、功耗散熱和電力限制、可管理性、編程簡易性、擴展性、管理維護費用等方面都面臨著巨大的挑戰(zhàn)。

1.2 PSTM研究現(xiàn)狀

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

1.3 GPU通用高性能編程技術

多核技術已經(jīng)代替單純CPU頻率的提升而成為計算機性能發(fā)展的主流。除了CPU技術之外，采用GPU、FPGA加速并與CPU異構并行的方式在這兩年大放異彩。使用這種架構的技術如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)的處理方法是把通用的計算問題轉換為GPU中能處理的圖形問題，然后通過圖形編程語言來進行編程。這樣的編程模式使得開發(fā)難度很高，開發(fā)者必須掌握強大的圖形處理能力。

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

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 分析結果

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

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

既然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等細粒度并行體系結構的加速。在此采用CUDA技術進行實現(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程序內核的執(zhí)行配置。根據(jù)GPU硬件資源，定義網(wǎng)格和線程塊，好的線程模型能大大提升程序的性能。

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

(2)線程塊（block）的定義：即定義一個線程塊有多少個線程和線程的組織方式，根據(jù)內核所需的寄存器數(shù)和共享內存數(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條線并行進行走時計算，可以使計算更均衡，性能更高。

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

圖2 異步IO邏輯結構圖

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

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

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

4 地震資料測試

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

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

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

4.2 性能結果

為了保證測試性能結果的穩(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運行結果的正確性。

通過對PSTM算法的性能剖析，分析算法的熱點函數(shù)，并對該熱點函數(shù)進行GPU并行化改造。初步改造移植結果說明，使用GPU進行加速可顯著提高疊前時間偏移計算速度，實驗數(shù)據(jù)證明了加速效果。測試數(shù)據(jù)表明，對于一個完整的PSTM體偏作業(yè)，一個GPU節(jié)點的運行時間是一個采用最新雙路CPU節(jié)點，并運行商用多線程CPU版PSTM時間的1/5。由此可見，PSTM高并行度、無數(shù)據(jù)依賴性等特征使得其較適合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].天然氣工業(yè)，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]李博.地震疊前時間偏移的一種圖形處理器提速實現(xiàn)方法[J].地球物理學報，2009，52(1)：245-252.