韓菲 李煒

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

關(guān)鍵詞：圖形處理器；疊前逆時(shí)偏移；混合粒度；并行計(jì)算；存儲(chǔ)優(yōu)化

中圖分類號(hào)：TP391 文獻(xiàn)標(biāo)識(shí)碼：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）

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

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

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

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

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

震源波場(chǎng)正演計(jì)算采用的聲波方程表示為式（1）

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

互相關(guān)成像條件式（4），由正演得到的震源波場(chǎng)S（x，y，z，t）和由逆時(shí)延拓得到的接收點(diǎn)波場(chǎng)R（x，y，z，t），利用互相關(guān)成像：

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

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

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

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

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

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

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

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

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

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

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

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

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

4.2.2 細(xì)粒度多線程并行

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

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

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

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

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

5.1 利用高速存儲(chǔ)器進(jìn)行優(yōu)化計(jì)算

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

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

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

5.2 利用常數(shù)存儲(chǔ)器進(jìn)行優(yōu)化

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

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

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

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

采用GPU集群進(jìn)行疊前逆時(shí)偏移并行運(yùn)算，一個(gè)節(jié)點(diǎn)上插有兩塊或四塊GPU卡，10個(gè)節(jié)點(diǎn)以上的GPU集群。節(jié)點(diǎn)之間使用MPI通信，節(jié)點(diǎn)內(nèi)部使用兩塊GPU進(jìn)行細(xì)粒度計(jì)算，保存的波場(chǎng)快照保存到當(dāng)前節(jié)點(diǎn)的磁盤上。

6 結(jié)論（Conclusion）

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

參考文獻(xiàn)（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.

作者簡(jiǎn)介：

韓 菲（1985-），女，博士，高級(jí)架構(gòu)師.研究領(lǐng)域：GPU并行計(jì)算，高性能計(jì)算，人工智能.

李 煒（1975-），男，碩士，高級(jí)架構(gòu)師.研究領(lǐng)域：高性能計(jì)算，高性能存儲(chǔ).