劉 帥 季曉慧* 蘆 俊 榮駿召

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

0 引言

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

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

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

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

1 多分量各向異性PSTM方法

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

(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波速度的各向異性參數(shù)。

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

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

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

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

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

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

2.1 CPU與GPU并行

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

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

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

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

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

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

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

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

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

圖2 并行算法流程

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

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

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

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

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

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

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

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

3 實際數(shù)據(jù)應(yīng)用與性能分析

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

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

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

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

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

圖5 縱波(a)與轉(zhuǎn)換波(b)偏移處理結(jié)果對比

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

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

表1 縱波與轉(zhuǎn)換波PSTM串、并行性能對比

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

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

4 結(jié)束語

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