基于MPI和CUDA的轉(zhuǎn)換波Kirchhoff疊前時(shí)間偏移并行計(jì)算

喻 勤，張少華，孔選林

（1.中國(guó)石油化工股份有限公司西南油氣分公司勘探開發(fā)研究院德陽分院，四川德陽 618000；2.中國(guó)石油化工股份有限公司多波地震技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川德陽 618000；3.東方地球物理公司科技信息處，河北涿州 072750）

轉(zhuǎn)換波Kirchhoff疊前時(shí)間偏移是轉(zhuǎn)換波地震數(shù)據(jù)處理中的一項(xiàng)關(guān)鍵技術(shù)，可實(shí)現(xiàn)全空間三維轉(zhuǎn)換波資料的準(zhǔn)確成像［1］。但轉(zhuǎn)換波勘探數(shù)據(jù)量巨大，進(jìn)行疊前偏移處理異常耗時(shí)。為了提高生產(chǎn)效率，目前大多數(shù)偏移算法采用大規(guī)模CPU 集群進(jìn)行計(jì)算，但因集群設(shè)備存在著成本高、功耗大、占用空間大以及維護(hù)成本高等缺點(diǎn)，其發(fā)展空間受到限制。近幾年，應(yīng)用GPU 高性能計(jì)算進(jìn)行地震偏移已成為一個(gè)新興的研究熱點(diǎn)，國(guó)內(nèi)外很多地球物理公司都開展了基于GPU 的偏移算法研究［2－5］，使疊前時(shí)間偏移、疊前深度偏移、單程波波動(dòng)方程偏移、逆時(shí)偏移等性能得到數(shù)十倍的提升。我們?cè)贚inux環(huán)境下，利用MPICH2.1，CUDA toolkit4.2，Eclipse等工具，實(shí)現(xiàn)了基于MPI和CUDA 的轉(zhuǎn)換波Kirchhoff疊前時(shí)間偏移并行算法。

1 轉(zhuǎn)換波Kirchhoff疊前時(shí)間偏移原理

在轉(zhuǎn)換波地震數(shù)據(jù)處理中，由于下行P 波和上行S波的路徑不對(duì)稱，需要做共轉(zhuǎn)換點(diǎn)（CCP）［6］選排，因此，轉(zhuǎn)換波資料處理流程一般為：抽CCP道集，速度分析，NMO，DMO，疊加和疊后偏移等。但因?yàn)镃CP 面元化不易確定轉(zhuǎn)換點(diǎn)真實(shí)位置，DMO 不能適應(yīng)層間速度劇烈變化和大陡傾角情況，所以需要進(jìn)行轉(zhuǎn)換波疊前時(shí)間偏移。該技術(shù)不需要抽取CCP道集和DMO 處理就能實(shí)現(xiàn)全空間三維轉(zhuǎn)換波資料的準(zhǔn)確成像。

根據(jù)各向異性雙平方根方程，疊前時(shí)間偏移方程［7－8］可以寫成

式中：tC為散射波旅行時(shí)；h是半源檢偏移距；x為炮檢距；ηeff，ζeff，γ0為各向異性參數(shù)；ΔtP是炮點(diǎn)到成像點(diǎn)的下行P波走時(shí)；ΔtS是成像點(diǎn)到檢波器的上行S波走時(shí)；vP2是下行P波速度；vS2是下行S波速度。由于從實(shí)際轉(zhuǎn)換波數(shù)據(jù)中很難獲得vS2，在各向異性速度分析時(shí)也得不到vP2，因此這兩個(gè)參量利用等效的轉(zhuǎn)換波速度vC2轉(zhuǎn)換得到。

圖1為轉(zhuǎn)換波成像示意圖，轉(zhuǎn)換波疊前時(shí)間偏移算法是計(jì)算出炮點(diǎn)到成像點(diǎn)的下行P 波走時(shí)tP，成像點(diǎn)到檢波器的上行S波走時(shí)tS，依據(jù)散射波旅行時(shí)公式（1），將上、下行波走時(shí)之和tPS上的能量累加到成像點(diǎn)處。轉(zhuǎn)換波Kirchhoff疊前時(shí)間偏移能將任何道集形式的輸入轉(zhuǎn)換成共成像點(diǎn)（CIP）道集。任意一個(gè)檢波點(diǎn)上接收的能量都可以依據(jù)炮檢點(diǎn)間的空間關(guān)系和各向異性旅行時(shí)方程分配到空間所有可能的成像位置，而所有炮檢對(duì)的能量都可以依據(jù)方程給出的射線路徑累加到圖1中散射點(diǎn)（成像點(diǎn)）的成像位置，從而使轉(zhuǎn)換波在三維空間任何位置準(zhǔn)確成像。

圖1 轉(zhuǎn)換波成像原理

轉(zhuǎn)換波疊前時(shí)間偏移從理論上規(guī)避了輸入數(shù)據(jù)零炮檢距的假設(shè)，從而避免了NMO 校正疊加產(chǎn)生的畸變，比疊后時(shí)間偏移能保存更多的疊前信息。疊前偏移后的疊加是共反射點(diǎn)的疊加，依據(jù)的模型是任意非水平層狀介質(zhì)，因此比疊后偏移成像在空間位置上更準(zhǔn)確。此外，疊前時(shí)間偏移采用均方根速度場(chǎng)，可以逐點(diǎn)進(jìn)行速度分析，成像精度較高，偏移速度較快，并可以提高AVO分析精度［9］。

2 基于MPI和CUDA 的算法實(shí)現(xiàn)

GPU 圖形處理器［10－13］是專門負(fù)責(zé)圖形渲染的核心處理器，內(nèi)部有上百個(gè)計(jì)算核心。CUDA 是Nvidia公司開發(fā)的一種SIMT 并行編程模型，用于開發(fā)高度并行的線程級(jí)程序。MPI［14－15］是一種基于消息傳遞并行編程模型的標(biāo)準(zhǔn)規(guī)范，從程序設(shè)計(jì)的角度理解是一個(gè)并行函數(shù)庫(kù)，其具體實(shí)現(xiàn)包括MPICH，OPEN MPI，LAM，IBM MPL 等多個(gè)版本，最常用和穩(wěn)定的是MPICH 和OPENMPI。MPI是目前超大規(guī)模并行計(jì)算應(yīng)用最廣、效率最高的并行編程模型，幾乎被所有并行環(huán)境（共享和分布式存儲(chǔ)并行機(jī)、集群系統(tǒng)等）和多進(jìn)程操作系統(tǒng)（UNIX，Linux，Windows NT）所支持，基于MPI開發(fā)的應(yīng)用程序具有可移植性、適用性、易用性、可擴(kuò)展性強(qiáng)和性能高等優(yōu)點(diǎn)。通過MPI可以方便地連接現(xiàn)有計(jì)算機(jī)，進(jìn)行高性能集群計(jì)算。

MPI與CUDA 的最大優(yōu)勢(shì)均在于并行計(jì)算，但是MPI面向的是主處理器CPU 進(jìn)程間的并行，屬于粗粒度級(jí)的并行，而CUDA 的并行是協(xié)處理器GPU 線程間的并行，屬于細(xì)粒度級(jí)的并行。要實(shí)現(xiàn)基于MPI和CUDA 的并行算法，需要對(duì)計(jì)算任務(wù)進(jìn)行劃分，偏向管理控制的部分采用進(jìn)程間的并行計(jì)算，計(jì)算密度大的部分采用線程間的并行計(jì)算。通過主處理器CPU 和協(xié)處理器GPU 結(jié)合，可以最大限度地發(fā)揮計(jì)算機(jī)的處理能力。

基于MPI和CUDA 的Kirchhoff疊前時(shí)間偏移并行算法設(shè)計(jì)主要包括以下幾個(gè)方面。

1）從硬件系統(tǒng)上考慮，采取一個(gè)CPU 綁定一個(gè)GPU 的方式，這樣可以形成最小計(jì)算粒度的節(jié)點(diǎn)，比較靈活地進(jìn)行數(shù)值計(jì)算，保證其單獨(dú)享有L2和L3緩存，避免多核CPU 產(chǎn)生的數(shù)據(jù)競(jìng)爭(zhēng)問題，減少資源需求，從而降低內(nèi)存負(fù)載，保持CPU 到顯存PCI－E信道、CPU 到主存以及主存到硬盤總線的高帶寬。

2）CPU 進(jìn)程間的通信通過MPI來完成，每個(gè)CPU 進(jìn)程綁定一個(gè)GPU。MPI主要用于完成均衡的任務(wù)分割，完成偏移過程所需相關(guān)參數(shù)的計(jì)算以及對(duì)各節(jié)點(diǎn)分發(fā)數(shù)據(jù)、管理控制各節(jié)點(diǎn)GPU 的計(jì)算流程等任務(wù)。每個(gè)節(jié)點(diǎn)的GPU 主要用于實(shí)際偏移計(jì)算，偏移計(jì)算所需的參數(shù)可以通過MPI分發(fā)到各個(gè)GPU 節(jié)點(diǎn)上。對(duì)速度數(shù)據(jù)建立了一個(gè)數(shù)據(jù)表，供各個(gè)GPU 快速查詢使用，通過疊加各個(gè)CPU 節(jié)點(diǎn)的計(jì)算結(jié)果，得出最終偏移剖面。偏移計(jì)算流程見圖2所示。

3）由于偏移本質(zhì)上是一個(gè)信號(hào)能量再分配的過程，即將接收信號(hào)的能量按照一定的關(guān)系分配到地層的各個(gè)反射點(diǎn)上，最終偏移后輸出的地震道和每一道檢波器所采集的數(shù)據(jù)都有關(guān)系，但相互之間是獨(dú)立的。因此可以把偏移過程簡(jiǎn)化為若干個(gè)獨(dú)立的、不相關(guān)的子問題，轉(zhuǎn)化為GPU 輕量級(jí)多線程應(yīng)用。

4）從GPU 硬件特性考慮，采取了多道輸入和多道輸出相結(jié)合的方式，一次讀取多道數(shù)據(jù)進(jìn)行偏移處理，得到整個(gè)成像空間的成像點(diǎn)數(shù)據(jù)再進(jìn)行輸出。該方式需要相對(duì)較少的輸入數(shù)據(jù)并加以重復(fù)利用，大幅度減少了主存和顯存之間的拷貝次數(shù)及拷貝時(shí)間，實(shí)現(xiàn)了數(shù)據(jù)的高吞吐并行計(jì)算。

5）設(shè)計(jì)中還需要考慮GPU 算法的均衡性，對(duì)于每個(gè)GPU 來說，偏移計(jì)算過程需要分塊（block）計(jì)算，每個(gè)塊中的計(jì)算需要線程同步，最為平均的計(jì)算量可以減少同步等待的時(shí)間。線程數(shù)量過多會(huì)導(dǎo)致線程同步時(shí)間變長(zhǎng)，也會(huì)導(dǎo)致計(jì)算過程中寄存器數(shù)量不足以及進(jìn)出堆棧帶來的損耗；線程數(shù)量過小，會(huì)導(dǎo)致塊數(shù)量增多，塊之間切換開銷相對(duì)較大。因此，在算法設(shè)計(jì)中，需要根據(jù)偏移計(jì)算需要的寄存器變量及計(jì)算過程的復(fù)雜性合理分配每個(gè)塊中的線程數(shù)量。

6）設(shè)計(jì)CUDA 偏移算法關(guān)鍵的一點(diǎn)，是如何在偏移計(jì)算過程中利用共用存儲(chǔ)器，在不發(fā)生存儲(chǔ)體沖突時(shí)，訪問延遲與寄存器相同。在不同的塊之間，共用存儲(chǔ)器是動(dòng)態(tài)分配的。在同一個(gè)塊內(nèi)，共用存儲(chǔ)器是進(jìn)行線程間低延遲數(shù)據(jù)通信的唯一手段。將偏移過程中每個(gè)塊內(nèi)的線程都需要的變量放入共用存儲(chǔ)器進(jìn)行訪問，例如最后的偏移能量累積疊加利用共用存儲(chǔ)器進(jìn)行，保證了訪問全局內(nèi)存的次數(shù)最少，避免了采用全局內(nèi)存的不合并訪問而造成的效率低問題，從而提高計(jì)算效率。

圖2 基于MPI和CUDA 的偏移計(jì)算流程

3 算法測(cè)試與效果分析

為測(cè)試基于MPI和CUDA 的Kirchhoff疊前時(shí)間偏移并行算法的性能，我們?cè)贚inux操作系統(tǒng)下實(shí)現(xiàn)了該算法。測(cè)試平臺(tái)為DELL 的6節(jié)點(diǎn)小型集群，集群中有2 個(gè)Tesla M2070 GPU 節(jié)點(diǎn)。測(cè)試環(huán)境配置如表1所示，測(cè)試數(shù)據(jù)為川西某區(qū)塊轉(zhuǎn)換波數(shù)據(jù)。選取其中4組進(jìn)行測(cè)試。

表1 測(cè)試環(huán)境

分別利用Intel Xeon CPU E5606單核，Tesla M2070GPU單卡以及MPICH2 和CUDA4.2 在集群上的2個(gè)GPU 節(jié)點(diǎn)對(duì)4組數(shù)據(jù)進(jìn)行了測(cè)試。以往采取的計(jì)算熱點(diǎn)部分比較忽略了如CPU 到GPU 數(shù)據(jù)拷貝傳輸過程、GPU 到CPU 的輸出過程以及其它一些會(huì)耗時(shí)的過程，雖然計(jì)算部分加速比非常高，但是整個(gè)計(jì)算任務(wù)的加速比低得多。為了獲得更客觀的計(jì)算效率，選取相同的計(jì)算量，即給出相同的計(jì)算任務(wù)，統(tǒng)計(jì)其以相同的偏移算法、不同的計(jì)算方式來完成所需要的總的時(shí)間（表2），最終獲得GPU 計(jì)算加速比，即串行執(zhí)行時(shí)間／并行執(zhí)行時(shí)間（表3）。這個(gè)加速比真正反映了一個(gè)計(jì)算任務(wù)的加速比。

從以上結(jié)果可以看出，基于CUDA 的偏移計(jì)算效率得到大幅度提升，即使在小數(shù)據(jù)量任務(wù)未能充分利用GPU 的情況下，計(jì)算效率也有大幅度提升。在數(shù)據(jù)量充分的情況下，GPU 單卡的速度比是CPU 單核的近200倍，利用MPI和GPU 的2個(gè)節(jié)點(diǎn)的計(jì)算速度比是CPU 單核的近400倍。針對(duì)4.76GB的測(cè)試數(shù)據(jù)，分別利用單核CPU，MPI和GPU 進(jìn)行偏移計(jì)算，得到的剖面如圖3 所示。將兩個(gè)剖面所對(duì)應(yīng)的點(diǎn)進(jìn)行數(shù)值誤差分析，結(jié)果如圖4所示，幅度誤差被控制在0.01的數(shù)量級(jí)，表明GPU 的計(jì)算精度對(duì)實(shí)際資料的成像結(jié)果基本沒有產(chǎn)生影響。

表2 測(cè)試計(jì)算時(shí)間 ms

表3 測(cè)試計(jì)算加速比

4 結(jié)束語

本文研究給出了基于MPI和CUDA 的轉(zhuǎn)換波Kirchhoff疊前時(shí)間偏移并行算法，測(cè)試分析結(jié)果表明：MPI和GPU（2 個(gè)節(jié)點(diǎn)）的計(jì)算速度是CPU（單核）的近400 倍，大幅度提高了轉(zhuǎn)換波Kirchhoff疊前時(shí)間偏移的計(jì)算效率，降低了計(jì)算成本。

基于MPI和GPU 的高性能計(jì)算技術(shù)有利于實(shí)現(xiàn)疊前時(shí)間偏移多次迭代，獲得合理匹配的速度模型，使得提高轉(zhuǎn)換波資料處理品質(zhì)成為可能。但是，GPU 偏移成像剖面還存在假頻現(xiàn)象，分辨率也有所降低，因此需要從反假頻和提高分辨率兩個(gè)方面對(duì)偏移算法進(jìn)一步進(jìn)行深入研究。

［1］馬昭軍，唐建明.疊前時(shí)間偏移在三維轉(zhuǎn)換波資料處理中的應(yīng)用［J］.石油物探，2007，46（2）：174－180 Ma Z J，Tang J M.The application of prestack time migration in 3Dconverted－wave data processing［J］.Geophysical Prospecting for Petroleum，2007，46（2）：174－180

［2］張兵，趙改善，黃俊，等.地震疊前深度偏移在CUDA 平臺(tái)上的實(shí)現(xiàn)［J］.勘探地球物理進(jìn)展，2008，31（6）：427－432 Zhang B，Zhao G S，Huang J.Implement prestack depth migration on CUDA platform［J］.Progress in Exploration Geophysics，2008，31（6）：427－432

［3］劉國(guó)峰，劉洪，王秀閩，等.Kirchhoff積分疊前時(shí)間偏移的兩種走時(shí)計(jì)算及并行算法［J］.地球物理學(xué)進(jìn)展，2009，24（1）：131－136 Liu G F，Liu H，Wang X M，et al.Two kinds of traveling time computation and parallel computing methods of Kirchhoff migration［J］.Progress in Geophysics，2009，24（1）：131－136

［4］李肯立，彭俊杰，周仕勇.基于CUDA 的Kirchhoff疊前時(shí)間偏移算法設(shè)計(jì)與實(shí)現(xiàn)［J］.計(jì)算機(jī)應(yīng)用研究，2009，26（12）：4474－4477 Li K L，Peng J J，Zhou S Y.Implement Kirchhoff prestack time migration algorithm on CUDA architecture［J］.Application Research of Computers，2009，26（12）：4474－4477

［5］劉國(guó)峰，劉欽，李博，等.油氣勘探地震資料處理GPU／CPU 協(xié)同 并行計(jì)算［J］.地 球 物理學(xué)進(jìn)展，2009，24（5）：1671－1678 Liu G F，Liu Q，Li B，et al.GPU／CPU co－processing parallel computation for seismic data processing in oil and gas exploration［J］.Progress in Geophysical（in Chinese），2009，24（5）：1671－1678

［6］Li X Y.Converted－wave moveout analysis revisited：the search for a standard approach［J］.Expanded Abstracts of 73rdAnnual Internat SEG Mtg，2003，805－808

［7］Li X Y，Dai H C，Mancini F.Converted－wave imaging in anisotropic media：theory and case studies［R］.Edinburgh：Edinburgh Anisotropy Project，2006

［8］Li X Y，Yuan J X.Converted wave traveltime equations in layered anisortropic media：an overview［R］.Edinburgh：EAP Annual Research Report，2001：3－32

［9］Dai H C，Li X Y，Conway P.3Dpre－stack Kichhoff time migration of PS－waves and migration velocity model building［J］.Expanded Abstracts of 74thAnnual Internat SEG Mtg，2004，1115－1118

［10］John D O.A surey of general－purpose computation on graphics hardware，eurographics ［EB／OL］.（2005－06－23）［2012－03－15］http：／／www.realworldtech.com／page.cfm？ArticleID＝RWT0908

［11］Nvidia.Nvidia CUDA compute unified device architecture programming guide version4.0［EB／OL］.［2011－12－9］http：／／developer.download.Nvidia.com／compute／cuda／4＿0／docs／CudaReferenceManual＿4.0.pdf

［12］張舒，褚艷麗，趙開勇，等.GPU 高性能計(jì)算之CUDA［M］.北京：中國(guó)水利電力出版社，2009：25－30 Zhang S，Chu Y L，Zhao K Y，et al.GPU high performance computing on CUDA［M］.Beijing：China Publishing Company of Water Conservancy and Electric Power，2009：25－30

［13］Budruk R，Anderson D，Shanley T.PCI express system architecture［EB／OL］.［2012－03－15］http：／／www.docin.com／touch／detail.do？id＝501507982

［14］都志輝.高性能計(jì)算之并行編程技術(shù)－MPI并行程序設(shè) 計(jì)［EB／OL］.（2007－10－9）［2012－03－15］http：／／wenku.baidu.com／view／00bacc204b35eefdc8d33375.html Du Z H.The parallel programme technolegy of high performance computing－MPI parallel programming［EB／OL］.（2007－10－9）［2012－03－16］http：／／wenku.baidu.com／view／00bacc204b35eefdc8d33375.html

［15］John D O，Mike H，David L，et al.GPU computing［J］.Proceedings of the IEEE，2008，96（5）：879－897