基于GPU的疊前逆時偏移混合粒度數(shù)據(jù)分割與存儲優(yōu)化

韓菲 李煒

摘 要：為提高疊前逆時偏移計算效率，本文采用MPI+CUDA混合粒度相結合的并行模式，對地震數(shù)據(jù)進行數(shù)據(jù)分割，合理劃分并行任務。總結出MPI+CUDA并行編程模型，提出疊前逆時偏移的混合粒度并行算法。根據(jù)CUDA特有的存儲方式，對疊前逆時偏移算法提出存儲優(yōu)化方案，更高效的利用GPU上各類存儲器，以進一步降低數(shù)據(jù)訪問所造成的時間延遲。

關鍵詞：圖形處理器；疊前逆時偏移；混合粒度；并行計算；存儲優(yōu)化

中圖分類號：TP391 文獻標識碼：A

Abstract：To improve the computational efficiency of prestack reverse-time migration，this paper adopts the MPI + CUDA parallel model to divide seismic data and parallel tasks.The MPI + CUDA parallel programming model is summarized and the hybrid granularity parallel algorithm of prestack reverse-time migration is proposed.Based on the special storage model of CUDA，we propose the storage optimization scheme for prestack reverse-time migration algorithm so as to reduce time delay caused by the data access with higher involvements of all kinds of memories on GPU.

Keywords：GPU；prestack RTM；hybrid granularity；parallel computing；storage optimization

1 引言（Introduction）

高性能計算技術在飛速發(fā)展，從微處理器由單核發(fā)展到多核，從并行機的出現(xiàn)到集群計算技術的飛速發(fā)展，無不證明計算機技術的飛速進步。當擁有多核處理器的GPU出現(xiàn)[1]，代替CPU處理大規(guī)模并行問題時，更是將高性能計算推向了一個新的高峰。目前，GPU已經廣泛應用于各個領域，并大大提高了計算密集問題的計算效率。

GPU在石油地球物理勘探領域有著廣泛應用。美國休斯頓的Headwave公司專門從事地學數(shù)據(jù)分析[2]，充分利用GPU的并行計算潛力，支撐實時工作流程、實時可視化和實時計算，支持疊前地震數(shù)據(jù)的分析與解釋。目前，GPU已在疊前逆時偏移中逐漸展開應用。因逆時偏移計算量且存儲量巨大，無法適應工業(yè)生產的需要，一直沒有在工業(yè)界得到廣泛應用，如何提高逆時偏移的計算效率成為逆時偏移能否大規(guī)模應用的重要環(huán)節(jié)[3]。如何合理地劃分計算任務，對計算數(shù)據(jù)進行合理存儲，提高存儲速度，成為我們首要研究的內容。本文采用粗、細粒度混合的數(shù)據(jù)分割方式劃分并行任務，給出疊前逆時偏移算法下MPI+CUDA的混合計算模型，研究了基于GPU疊前逆時偏移的并行算法。根據(jù)CUDA特有的存儲方式，對逆時偏移算法的存儲提出優(yōu)化方案。

2 疊前逆時偏移原理（The principle of prestack reverse time migration）

疊前逆時深度偏移求解兩次波動方程，即從源點以時間正方向正演出空間所有點的波場及從檢波點位置以時間逆序推出所有空間反射點的波場；利用相應成像條件對兩個波場進行成像。逆時偏移是全波場的雙程波動方程偏移方法[4]。可以對復雜速度體成像，不受成像傾角限制（可達到180°）；能對所有波動現(xiàn)象（如反射波、折射波、繞射波、回轉波、棱柱波和多次波等）進行成像。此過程中波的傳播是相對完整的，因而所用偏移速度模型無論如何復雜都不需要平滑處理，但對于常規(guī)Kirchhoff偏移是很難做到的。

首先，給出該方法的具體實施流程：（1）以時間正方向正演計算震源波場；（2）以時間逆方向逆推接收點波場；（3）利用成像條件，對兩波場成像；（4）所有炮集成像結果疊加得到最后成像結果。

震源波場正演計算采用的聲波方程表示為式（1）

對于式（2）和式（3）的數(shù)值求解方法較多，如有限差法分、有限元和偽譜法等，每種方法都有各自的優(yōu)缺點。本文采用的數(shù)值求解方法是時間方向為四階差分精度的偽譜法。通過數(shù)值計算得到震源波場pF和接收點波場pB后，然后利用相應成像條件對兩個波場進行成像。那么，不同成像條件的應用，將直接影響成像的振幅、分辨率及其保幅性。所以，成像條件的選擇比較關鍵，對保幅性的疊前深度偏移更是至關重要。

互相關成像條件式（4），由正演得到的震源波場S（x，y，z，t）和由逆時延拓得到的接收點波場R（x，y，z，t），利用互相關成像：

3 CPU/GPU協(xié)同架構（CPU/GPU collaborative architecture）

CUDA是一種將GPU作為數(shù)據(jù)并行計算設備的軟硬件體系，目前已廣泛應用。在運算過程中GPU（設備）是作為CPU（主機）的協(xié)處理器來執(zhí)行操作的，主機和設備均有自己的存儲器。CPU-GPU并行編程模型如圖1所示[5]。

4 混合粒度數(shù)據(jù)并行算法（Hybrid-grained data parallel algorithm）

粒度是對一個并行計算任務中計算量大小和通信量大小比值的度量。將問題分解為多個子問題，根據(jù)并行計算任務的大小，可以分為粗粒度并行、細粒度并行和中粒度并行。本文采取粗粒度并行與細粒度并行兩種并行方式設計算法。

4.1 粗粒度數(shù)據(jù)并行

4.1.1 粗粒度數(shù)據(jù)分割

粗粒度是將一個任務分解為含有較多計算量的多個子問題，這些子問題可以獨立地并行解決問題。本文采用MPI粗粒度并行方式對地震數(shù)據(jù)進行數(shù)據(jù)分割。首先由主節(jié)點（Master node）根據(jù)有效節(jié)點數(shù)Nn和地震數(shù)據(jù)總炮數(shù)Ns進行節(jié)點間作業(yè)任務分配。為了提高讀取數(shù)據(jù)的速度，采用多節(jié)點并行讀取數(shù)據(jù)，能大大提高數(shù)據(jù)輸入輸出效率。每個節(jié)點單炮偏移完后，返回主節(jié)點判斷是否所有炮集數(shù)據(jù)都計算處理完畢，若沒完成進入下一次循環(huán)，直至所有炮集資料都偏移完為止。粗粒度數(shù)據(jù)分割如圖2所示。

4.1.2 粗粒度數(shù)據(jù)并行

利用多GPUs和多核CPU協(xié)同加速計算逆時偏移，主要思想：因逆時偏移是按照炮集一炮一炮進行偏移，采用MPI啟動多進程，每個進程讀取一炮地震數(shù)據(jù)，將一炮地震數(shù)據(jù)發(fā)送到一塊GPU上進行運算，在GPU上啟動多線程并行執(zhí)行一炮地震數(shù)據(jù)的偏移過程。通過調度實現(xiàn)在多GPUs上并行，將讀取文件和其他步驟設計成串行過程。充分利用GPU高效的浮點數(shù)處理能力和多核CPU良好的任務分配和調度能力。多GPUs算法流程如圖3所示。

4.2 細粒度數(shù)據(jù)并行

4.2.1 細粒度數(shù)據(jù)分割

細粒度是一個任務劃分成包含較短的程序段和較小的計算量的并行任務，粒度越小，越可開發(fā)更多的并行性，提高并行度，這是有利的方面。但粒度越小，通信次數(shù)和通信量就相對增多，增加了額外開銷。GPU可以同時開啟成百上千個線程，擁有強大的數(shù)據(jù)處理能力，專用于解決數(shù)據(jù)并行計算的問題，具有極高的計算密度（數(shù)學運算與存儲器運算的比率），是典型的細粒度并行編程模型。

本文采用GPU上細粒度數(shù)據(jù)分割方式，將MPI一個進程讀取的一炮地震數(shù)據(jù)通過主機發(fā)送到GPU上。每炮地震數(shù)據(jù)有104道地震記錄，每道地震記錄發(fā)送到GPU的一個線程上，每個線程計算一道地震數(shù)據(jù)。線程細粒度數(shù)據(jù)分割如圖4所示。

4.2.2 細粒度多線程并行

在GPU上啟動多線程，將一炮地震數(shù)據(jù)中104道地震數(shù)據(jù)映射到線程中，以取代正演和逆推過程中偏移一炮地震數(shù)據(jù)從第一道循環(huán)到第104道，從第1個采樣點到1024采樣點的兩層循環(huán)。將計算量均勻地加載到各線程，提高算法的并行度，使算術邏輯單元盡可能均勻地分配到計算任務，是提高GPU計算效率的舉措之一。

4.3 疊前逆時偏移MPI+CUDA混合粒度并行算法

結合粗粒度與細粒度并行的優(yōu)點，考慮采用MPI+CUDA多粒度混合編程模型，節(jié)點間通過MPI消息傳遞機制，實現(xiàn)粗粒度并行；節(jié)點內利用GPU強大的并行數(shù)據(jù)處理能力，使用CUDA架構啟用多線程機制，利用共享內存實現(xiàn)高效快速細粒度的數(shù)據(jù)并行和線程并行?；旌狭６染幊棠Ｐ腿鐖D5所示。MPI+CUDA混合編程模型能充分發(fā)揮集群節(jié)點間分布式存儲和節(jié)點內共享存儲的優(yōu)勢，既可以發(fā)揮每個節(jié)點的巨大計算能力，又可以充分利用集群的可擴展性，使并行效率顯著提高，為CPU+GPU混合架構做大型計算提供了一種有效的并行策略。

疊前逆時偏移算法的主要部分：震源波場的正傳，檢波點波場的反傳，以及相關成像條件均是CPU串行計算最為耗時的部分。針對共炮域的逆時偏移而言，本文采取混合粒度的并行方式，以單個炮集為粗粒度并行，通過GPU并行計算的方法，將每個炮集包含的104道地震數(shù)據(jù)發(fā)送到GPU中執(zhí)行線程級細粒度并行。在GPU上，先將速度數(shù)據(jù)、一炮的地震數(shù)據(jù)和激發(fā)點的初始波場值由內部存儲器（內存）傳至設備存儲器（顯存），這樣可以避免由多次重復對內部存儲器的數(shù)據(jù)讀寫所引起的時間延遲。然后把GPU的多核處理器劃分為相應個數(shù)的計算塊（block），同時每個計算塊又可以劃分為若干個線程（thread），利用GPU的多線程實現(xiàn)大規(guī)模的并行計算。這樣就可以使得在計算過程中涉及的數(shù)據(jù)讀寫和運算均在設備存儲器和圖形處理器（GPU）中進行，將波場延拓及相關成像的計算并行化，達到提高計算效率的目的。最后再將數(shù)據(jù)傳回內部存儲器，通過CPU進行結果的I/O操作。

5 逆時偏移存儲優(yōu)化（Reverse time migration storage optimization）

5.1 利用高速存儲器進行優(yōu)化計算

寄存器是多處理器上最快的片上存儲器。疊前逆時偏移并行算法內核函數(shù)中的變量在允許情況下都定義為寄存器型。每個SM具有32位寄存器，如果內核編譯時使用了超過執(zhí)行配置允許數(shù)量的寄存器，會造成內核由于寄存器不足而無法啟動。疊前逆時偏移并行算法中每個SM使用寄存器小于32kB，沒有超限。

CUDA提供了具有高速讀寫訪問并行數(shù)據(jù)的共享存儲器，其速度在無存儲體沖突時等同于寄存器。共享存儲器是可以被同一塊中的所有線程訪問的可讀寫存儲器，應用程序可以利用它來最小化對DRAM的過度提取和巡回，從而降低對DRAM存儲器帶寬的依賴程度。并行算法內核函數(shù)中除了部分計算參數(shù)和中間型數(shù)據(jù)運算時需要共享型外，把用于接收從設備存儲器傳入的目標區(qū)和搜索區(qū)的地震數(shù)據(jù)也定義為共享型。由于每個SM上只有48kB共享存儲器，一個SM可以同時執(zhí)行多個線程塊，片內的存儲資源則會被分配給這些并發(fā)的線程塊中。根據(jù)每個SM的資源上限，計算分發(fā)單元就可以知道最多能在一個SM上分配多少個block。若在計算相關系數(shù)的內核函數(shù)kernel中每個block中有104個線程，每個block使用49152B共享存儲器，每個線程使用四個寄存器，那么可以有：每個block使用的共享存儲器：49152Byte；每個block使用的寄存器數(shù)量：4*104=416；每個block中的warp數(shù)量：ceil（416/32）=13。

根據(jù)每個block使用的資源，就可以計算出kernel中由每個因素限制的最大活動塊數(shù)量：每個SM中的最大活動塊數(shù)量：8；由共享存儲器數(shù)量限制的活動塊數(shù)量：floor（393216/49152）=8；由warp數(shù)量限制的活動塊數(shù)量：floor（32/13）=2。

5.2 利用常數(shù)存儲器進行優(yōu)化

對當前硬件來說，一個half-warp中的線程訪問常數(shù)存儲器中數(shù)據(jù)，僅需要1—100個時鐘周期就可以獲得這個數(shù)據(jù)。實際使用常數(shù)存儲器時的速度一般還是低于寄存器或者共享存儲器，但還是明顯高于將數(shù)據(jù)存放在全局存儲器中的情況。常數(shù)存儲器空間較小，只有64kB。由于疊前逆時偏移在算法中還需讀取炮點坐標文件和每一炮含有104道的計數(shù)文件，大小為2.9kB和960Bytes，而且是相對不變量，因此，在疊前逆時偏移并行算法中可將炮點坐標文件和統(tǒng)計每炮含有多少地震道文件放入常數(shù)存儲器中，用來加速坐標的讀取。

5.3 全局存儲器存儲優(yōu)化

疊前逆時偏移算法中在震源波場順時外推的過程中需要保存波場快照，若以每5ms保存一次快照來計算，一炮數(shù)據(jù)偏移過程中保存的波場快照約為3.5G，原有的算法將如此大的數(shù)據(jù)量保存到刀片機當前節(jié)點的硬盤上，不同進程將不同數(shù)據(jù)保存在不同的節(jié)點硬盤上，同時進行多節(jié)點I/O存取，時間延遲會非常大。本文選用Tesla C2075 GPU型號，該GPU的全局存儲器的容量達到4GB，因此可將單炮偏移結果保存到GPU的全局存儲器下，然后直接進行檢波點波場逆推，省去原來從磁盤讀到內存，從內存讀到CPU的I/O傳輸時間，同時也省去CPU-GPU主機與設備之間低吞吐量數(shù)據(jù)傳輸?shù)臅r間。

5.4 GPU集群存儲優(yōu)化

采用GPU集群進行疊前逆時偏移并行運算，一個節(jié)點上插有兩塊或四塊GPU卡，10個節(jié)點以上的GPU集群。節(jié)點之間使用MPI通信，節(jié)點內部使用兩塊GPU進行細粒度計算，保存的波場快照保存到當前節(jié)點的磁盤上。

6 結論（Conclusion）

為提高疊前逆時偏移計算效率，本文運用MPI+CUDA混合粒度相結合的并行方式，提出疊前逆時偏移算法中更加合理的地震數(shù)據(jù)分割方式，更加合理地劃分并行任務，并總結出MPI+CUDA并行編程模型，將此模型應用于疊前逆時偏移算法中，提出疊前逆時偏移的混合粒度并行算法。同時，本文還提出疊前逆時偏移算法的存儲優(yōu)化策略，更高效的利用GPU上各類存儲器，以進一步降低數(shù)據(jù)訪問所造成的時間延遲。

參考文獻（References）

[1] Foley D，Danskin J.Ultra-Performance Pascal GPU and NVLink Interconnect[J].IEEE Micro，2017，37（2）：7-17.

[2] Marchesini S，Krishnan H，Daurer B J，et al.SHARP：a distributed GPU-based ptychographic solver[J].Journal of Applied Crystallography，2016，49（4）：1245-1252.

[3] Tan Y，Ding K.A Survey on GPU-Based Implementation of Swarm Intelligence Algorithms[J].IEEE Transactions on Cybernetics，2016，46（9）：2028-2041.

[4] Obrecht C，Asinari P，Kuznik F，et al.Thermal link-wise artificial compressibility method：GPU implementation and validation of a double-population model[J].Computers & Mathematics with Applications，2016，72（2）：375-385.

[5] Wang T，Jiang Z，Kemao Q，et al.GPU Accelerated Digital Volume Correlation[J].Experimental Mechanics，2016，56（2）：1-13.

作者簡介：

韓 菲（1985-），女，博士，高級架構師.研究領域：GPU并行計算，高性能計算，人工智能.

李 煒（1975-），男，碩士，高級架構師.研究領域：高性能計算，高性能存儲.