CPU與GPU協(xié)同并行的多分量地震數據各向異性疊前時間偏移

劉 帥 季曉慧* 蘆 俊 榮駿召

(①中國地質大學(北京)信息工程學院，北京 100083；②中國地質大學(北京)能源學院，北京 100083)

0 引言

傳統(tǒng)地震勘探利用縱波信息解決構造問題，但隨著剩余油氣藏勘探難度的增加，僅使用縱波信息已不能滿足現代油氣勘探的需求。多分量地震勘探技術可以充分利用縱波和轉換波信息[1-2]。

隨著各向異性理論的發(fā)展，多分量地震數據的各向異性疊前時間偏移(Pre-stack time migration，PSTM)已成為地震勘探的前沿技術，各向異性參數也被證明是提高地震偏移成像精度的關鍵因素[3]。Kirchhoff疊前時間偏移算法是經典偏移算法，該算法計算量龐大，也是整個常規(guī)數據處理周期中最為費時的環(huán)節(jié)[4-5]。隨著勘探目標數據體越來越大，地震數據處理周期也越來越長，難以滿足地震資料高效解釋的需求，亟待嘗試并行計算縮短處理時間，提高計算效率[6-7]。

前人已經進行了一些關于疊前時間偏移算法的并行化研究，包括基于多節(jié)點CPU的MPI(Message passing interface，消息傳遞接口)技術和基于GPU(Graphics processing unit，圖形處理器)的CUDA(Compute unified device architecture，計算同一設備架構)并行架構。其中，馬召貴等[8]提出了一種基于CPU的多級并行優(yōu)化方案，并應用于起伏地表的疊前時間偏移中。該方案消除了成像域并行的通訊耗費，降低了系統(tǒng)I/O對計算性能的影響。Alves等[9]提出一種針對CPU指令集的手動向量化算法優(yōu)化方法，在不使用協(xié)處理的情況下，算法加速比達到了4倍。但是CPU限于其自身的硬件結構特點，計算密集型算法仍不適宜在CPU結構上展開。李肯立等[10]利用CUDA將GPU協(xié)同CPU計算應用到縱波Kirchhoff疊前時間偏移，使用單個8800GT型號GPU協(xié)同CPU并行，較CPU串行算法獲得了6倍的加速效果，但其使用測試數據的規(guī)模及計算效果很難滿足實際應用的需求。喻勤等[11]提出基于CUDA和MPI的多節(jié)點并行Kirchhoff疊前時間偏移算法，使用兩個M2070型號GPU節(jié)點并行，加速比達到了300。使用多個節(jié)點必然會增加服務器成本，并且導致通信開銷增大[12]，如果通過增加GPU個數增強單個服務器的計算能力，那么在節(jié)省部分成本的同時，還可達到與多節(jié)點類似的加速效果，并具有較低能耗[13]。

本文使用縱波、轉換波數據進行疊前時間偏移，全面地反映地下復雜構造，并且加入各向異性參數提高成像精度；同時將OpenMP和CUDA結合在一起，實現單節(jié)點上CPU與多個GPU協(xié)同并行的疊前時間偏移算法，使用內存映射的方法進行地震數據的I/O優(yōu)化以及多主測線(Inline)輸入輸出的方法讀寫數據，在降低計算成本的同時，使大規(guī)模地震數據資料可以在小內存服務器上并行計算，達到了較高加速比。

1 多分量各向異性PSTM方法

本文算法是基于Kirchhoff積分各向異性疊前時間偏移算法，其流程如圖1所示，即將特定位置的輸入地震數據(PP波為CMP道集，PS波為ACP道集)疊加映射到成像空間(CIP道集)中。其中，數據位置由計算地震波上、下行旅行時得出。相比常規(guī)疊前時間偏移，各向異性疊前時間偏移在每個網格點上加入了反映當前工區(qū)點地下特征的各向異性參數η，以提高成像精度，因此計算量會有常量級的增加。具體旅行時計算方法由楊震等[14]根據文獻[15]改進得到。即有VTI介質精確時差公式

(1)

(2)

式中：tPS是PS波旅行時；tPP是PP波旅行時；xP和xP1分別是PS波和PP波的下行P波的炮點到成像點水平距離；xS和xP2分別是PS波上行S波和PP波上行P波的檢波點到成像點水平距離；t0PP和t0PS分別是PP波與PS波自激自收時間；VP是PP波NMO速度；VS是S波NMO速度；θP和θP1分別是PS波和PP波入射角；θS和θP2分別是PS、PP波的出射角；ηP和ηS分別是P波和S波速度的各向異性參數。

如圖1所示，串行算法采用面向成像域的實現策略，即先讀取工區(qū)屬性和偏移參數，根據讀到的成像主測線、聯(lián)絡線范圍及測線間隔建立成像網格。在成像網格中按順序選定一道成像道(一道成像道即一個成像網格點)，根據成像道位置讀取一道速度數據，根據成像道位置及偏移參數中的擴展線數讀取一道輸入地震數據，其中，擴展線數x作用是當偏移第A(數字)號線時，輸入A-x到A+x號線地震數據共同對第A號線進行偏移，以提高偏移成像效果及分辨率。

圖1 疊前時間偏移串行算法流程

同時，算法的整體計算量相比無擴展線數算法成倍增加。從成像道內按順序選定一個采樣點，將該點位置、偏移參數及讀取的速度數據等代入式(1)或式(2)后，計算該道輸入地震數據相對該采樣點的地震波旅行時。 然后提取該道輸入地震數據中，計算所得旅行時對應位置的地震數據，并疊加到該采樣點上，第一輪判斷成像道內每個采樣點計算是否完畢：“否”則選取成像道內下一個采樣點進行偏移計算。“是”則第二輪判斷是否將擴展線內全部輸入地震數據讀取完畢：“否”則讀取擴展線內下一道輸入地震數據，“是”則第三輪判斷是否全部成像道計算完畢：“否”則選取成像網格中下一道成像道，“是”則輸出成像空間，算法結束。

從圖1可知，算法的計算核心是計算地震波旅行時及輸入數據疊加采樣點部分，疊前時間偏移算法整體的計算量取決于成像網格點數k、擴展線中包含的輸入地震道數m及每道數據的采樣點數n，算法時間復雜度為O(kmn)。如果對由上百萬(106)地震道組成、每道采樣點數為103的小規(guī)模數據以擴展線數為2(擴展線包含的輸入地震道數約為104)、成像網格點數為104的成像參數進行偏移，算法的時間復雜度約為1011，需在服務器上運行數十小時。需要對其進行并行化運算，而算法并行化的必要條件是計算過程中不存在數據依賴，即成像網格內任一成像點的偏移不依賴于其他成像點。從圖1中可以看出，本算法中的各個成像點之間相互獨立，不存在數據依賴，具備并行化條件。

除計算量巨大外，運算過程內存消耗也很大。實際地震數據處理時，無法將地震數據及速度數據整體存放于內存中，因此必須調整計算思路。

2 CPU與GPU協(xié)同并行偏移算法

2.1 CPU與GPU并行

并行指在同一時間內由多個處理器同時處理同一個任務的不同部分，從而達到加速處理的效果[16-17]。

CPU和GPU都可以通過成熟的編程框架進行并行編程，但CPU的結構設計偏向于計算機邏輯指令的編譯和執(zhí)行，而GPU的結構設計更偏向于大規(guī)模的數據計算，其浮點數運算速度可達到CPU的數十倍[18-20]，因此合理使用并行框架，發(fā)揮各自擅長的能力處理問題成為并行化的關鍵。OpenMP和CUDA分別是常用基于CPU和GPU的并行框架。

OpenMP是一種用于共享內存系統(tǒng)的CPU多線程并行方案，提供了很多庫函數來解決常見的多線程程序設計問題[21-22]。在本文中，其作用是在CPU上開啟多個線程與多個GPU一一對應，將讀取到的數據分發(fā)給不同的GPU并控制對應的GPU進行數據計算，以及計算完成后接收多個GPU的計算結果，按順序寫入結果文件。

CUDA是NVIDIA(英偉達)公司發(fā)布的應用于GPU的并行編程模型。在CUDA中，GPU及其內存稱為設備(Device)，CPU以及系統(tǒng)的內存稱為主機(Host)。由于主機內存與設備內存相互獨立，所以，在實際使用CUDA時，必須將GPU計算所需數據從主機內存?zhèn)鬏數皆O備內存，以便GPU讀取數據進行計算，且在計算完成后，必須將計算結果從設備內存?zhèn)鬏數街鳈C內存，以便CPU輸出計算結果。使用CUDA提供的類C語言函數就可實現分配設備空間、主機和設備間的數據傳輸以及釋放設備空間等必要操作。

設備內存包括幾種不同特點的存儲結構，根據不同數據的讀取方式，使用與之適應的存儲方法，以提高并行算法性能。全局內存(Global memory)是可讀寫的設備內存，容量最大，但是讀取速度慢，用來存儲從主機傳輸到設備上的數據。紋理內存(Texture memory)和常量內存(Constant memory)是在核函數執(zhí)行期間只讀的設備內存，分別緩存在紋理緩存(Texture cache)和常量緩存(Constant cache)上，容量小但是讀取速度快。此外，紋理內存的特點是當多個線程在讀取數據的位置非常鄰近時，可以加快數據的讀取速度。

核(Kernel)函數是每個GPU線程在GPU設備上執(zhí)行的函數，一般對應串行算法中計算最密集、計算量最大的部分。CUDA通過控制GPU中的流處理器(Stream processor，SP)，同時啟動大量GPU線程，執(zhí)行同一個核函數對數據進行操作，達到并行計算的目的。在啟動CUDA核函數進行計算時，會伴隨著分配設備存儲空間、主機與設備間相互傳輸數據以及回收設備上的存儲空間等一系列必要操作以配合計算，這些操作存在一定的時間開銷，所以在設計并行算法時，應該盡量減少核函數啟動次數，降低花費在分配空間、傳輸數據以及回收空間操作上的時間[23-24]。

3.2 CPU與GPU協(xié)同并行偏移算法實現

本文提出的CPU與多個GPU協(xié)同并行算法的流程圖如圖2所示。

本算法將全部偏移任務分成多個計算周期完成，點線框包含的步驟為一個計算周期，由CPU進行算法中的邏輯處理、數據讀寫及分發(fā)，由各GPU進行旅行時計算，并根據得到的旅行時偏移輸入地震數據。

圖2 并行算法流程

具體方案是由每個GPU負責一條主測線的成像，因此一個計算周期內可以偏移計算GPU個數條成像主測線。在GPU內部，每個GPU線程與讀取的擴展線內每道輸入地震數據一一對應，以控制各道輸入地震數據的偏移。各GPU線程并行計算各道輸入地震數據相對各成像道(一個成像網格點)內每個采樣點的地震波旅行時，再根據計算所得旅行時，從輸入地震數據中提取數據，疊加到相應的成像采樣點上。本文采用在擴展線內包含的輸入地震數據道數粒度(104)并行的原因，主要是該粒度數量級與GPU中可同時并發(fā)線程數量級一致，每個GPU線程的計算量比較均衡，能充分利用GPU多線程并發(fā)的特點及GPU的浮點計算能力。

本文以包含N條主測線(線號為1～N)的地震數據，使用4個GPU協(xié)同CPU并行成像計算為例，解釋算法流程；同時，考慮到計算效率和成像效果，將擴展線數設置為2。

算法的第一個計算周期在GPU上的并行計算方法如圖3所示，首先將計算所需的地震數據與速度數據由硬盤讀入CPU內存，在第一個計算周期，即對1到4號測線進行偏移計算時，由于成像空間位于工區(qū)邊界，根據擴展線數為2，需讀取1到6號共6條主測線數據及速度模型數據。由于地震數據及速度數據量較大，本文在CPU端使用內存映射方法進行讀取，即把地震數據文件映射進CPU內存，得到一個指向整體數據起始位置的指針，通過給指針附加偏移量讀取指定位置的數據。并建立一個地址索引數組來記錄各地震道相對起始位置的偏移量，以幫助算法快速而準確地切換數據讀取位置。速度模型數據使用相同的內存映射方法進行讀取。

CPU讀取相關數據后，利用OpenMP開啟多個CPU線程，每個CPU線程控制一個GPU設備，各個CPU線程將讀取到的數據分塊傳輸到各個GPU上，過程如圖3的中部所示。由于第一個周期的測線處在成像空間左邊界且擴展線數為2，因此只有1到3號主測線數據傳輸到0號GPU上偏移計算1號成像域主測線，1到4號主測線數據傳輸到1號GPU上偏移計算2號成像域主測線。其他非工區(qū)邊界情況如圖3中2號GPU與3號GPU部分所示。因擴展線數為2，故使用5條輸入主測線數據偏移計算1條成像主測線。對于量大的地震數據，將其傳輸到空間較大但讀取相對較慢的GPU全局內存上，對于頻繁讀取且讀取位置鄰近的速度數據，本文將其綁定在讀取較快但空間相對較小的GPU紋理內存上。

各GPU接收到數據后，如前所述每個GPU負責一條成像域主測線的偏移成像計算，如圖3所示，0號GPU偏移計算1號成像域主測線的成像空間，1號GPU偏移計算2號成像域主測線的成像空間，以此類推。然后以一個GPU線程對應一道地震道數據的原則，開啟各GPU線程并啟動GPU上的核函數并行偏移計算。通過每個線程專屬的線程號控制核函數，依次計算各地震道在不同炮檢距的地震波的旅行時，然后各GPU線程將旅行時對應的輸入地震數據疊加到成像采樣點上。各GPU內成像主測線并行偏移計算完畢后，各GPU分別傳輸各自所得成像空間數據到CPU，并且釋放GPU上所有分配的內存空間。所得成像空間數據由CPU按順序寫入到偏移結果文件中，本計算周期計算結束。

第二個計算周期與第一個計算周期過程相似，不同點在于第二個計算周期將針對5～8號成像空間進行并行偏移計算，由于不涉及工區(qū)邊界，所以需讀取3～10號共8條輸入測線數據和速度數據，分發(fā)數據時也按照各GPU偏移計算所需數據范圍進行分發(fā)，依次循環(huán)，直至全部成像空間計算完畢。

圖3 并行算法第一個計算周期示意圖

當運算到工區(qū)右邊界時，若剩余成像空間主測線數小于當前開啟GPU個數，則自動修正參與計算的GPU個數為剩余成像空間主測線數，以保證一個GPU偏移計算一條成像域主測線。

3 實際數據應用與性能分析

測試數據來自新場工區(qū)的X分量和Z分量數據，主測線范圍是1380～1595，聯(lián)絡線范圍是1002～1779，每個地震道采樣點數為3500，采樣點間隔為2ms，每個分量數據量約為29GB。成像網格主測線范圍是1380～1595，聯(lián)絡線范圍是1002～1779，成像的時間長度為6998。硬件環(huán)境方面，CPU為2路Intel Xeon E5-2695v2，核心頻率為2.4GHz，設計功耗為115W。每個GPU為NVIDIA Tesla K40M計算卡，該型號有2880個CUDA計算核心，存儲空間約為12GB，存儲帶寬為288GB/s，單精度浮點數計算峰值為4.29Tflops(每秒萬億次浮點運算)，雙精度浮點數計算峰值為1.43Tflops，設計功耗為235W。

首先對比CPU串行算法與本并行算法結果，以轉換波疊前時間偏移算法為例，隨機提取某一相同位置剖面(圖4)，左側為CPU串行計算結果，右側為并行算法計算結果，可見反映的地下構造幾乎相同，證明了并行算法的正確性。

圖5分別是縱波與轉換波使用本文并行PSTM算法偏移結果的同一測線剖面，從波組特征可見縱波2.5～3.1s層位(紅箭頭)對應轉換波3.5～4.5s時段，縱波3.6～3.9s層位(藍箭頭)對應轉換波5.0～5.5s時段。

表1是PSTM串行算法與不同GPU個數協(xié)同CPU算法并行的運行時耗統(tǒng)計，從加速比可看出并行算法的性能指標明顯較高。圖6是縱波PSTM算法的并行效果與理想加速效果的折線圖對比，其中，理想加速比用來判斷CPU的并行效果，是指在使用OpenMP開啟多CPU線程控制多GPU協(xié)同CPU計算時，各CPU線程之間無線程同步等損耗的線性加速比，即實際加速比曲線越接近線性理想加速比，算法的并行度越高，并行效果越好。

圖4 CPU串行(a)、并行(b)算法偏移處理結果對比

圖5 縱波(a)與轉換波(b)偏移處理結果對比

可見通過使用OpenMP技術，同時使用多個GPU偏移計算，具有較好加速效果，未達到線性加速比的原因主要是由于隨著CPU開啟線程數的增加，在線程同步上花費的時間也隨之增加。由于轉換波PSTM相比縱波PSTM計算量稍多，所以實際加速效果也稍好。

內存使用方面，由于本算法分割了成像空間，所以在每個計算周期，系統(tǒng)的內存只存放該計算周期所需數據，在原始地震數據成倍增長的情況下，也不會導致內存溢出、無法計算的情形。

表1 縱波與轉換波PSTM串、并行性能對比

圖6 縱波(a)及轉換波(b)PSTM并行方式實際加速比與理想加速比對比

此外，現在業(yè)界普遍使用CPU集群加速計算，而CPU并行試算的理想加速比是并行使用的核心總數，但由于并行模型的通信等問題，實際加速比往往低于理論加速比，并且本文提出的6個GPU協(xié)同CPU的并行算法加速比已能達到444和449，若想通過CPU集群計算達到類似加速比，則至少需要450個CPU核心的計算資源，相當于一個數+計算節(jié)點的CPU集群規(guī)模。所以從系統(tǒng)造價及能耗角度看，GPU并行明顯優(yōu)于CPU并行。

4 結束語

本文提出了基于CPU與GPU協(xié)同并行的多分量地震數據各向異性疊前時間偏移算法，優(yōu)化了傳統(tǒng)偏移方法的I/O方式，使用多次讀取、每次讀取多主測線、每次計算GPU個數條成像主測線的方法分解計算任務，兼顧了算法占用的內存空間和計算效率。在一臺服務器上對實際工區(qū)約29G的縱波及轉換波數據使用6個GPU協(xié)同CPU進行PSTM算法并行時，縱波PSTM算法的加速比可以達到444，轉換波PSTM算法的加速比可以達到449，加速比相比傳統(tǒng)的通過增加CPU集群數量的并行方式有大幅度提高，且充分利用了單個服務器上的設備接口數量，相比多個單GPU計算節(jié)點，避免了節(jié)點間的通信開銷同時降低了計算成本。