喻 勤，張少華，孔選林

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

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

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

在轉(zhuǎn)換波地震數(shù)據(jù)處理中，由于下行P 波和上行S波的路徑不對稱，需要做共轉(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ù)不需要抽取CCP道集和DMO 處理就能實(shí)現(xiàn)全空間三維轉(zhuǎn)換波資料的準(zhǔn)確成像。

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

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

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

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

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

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

GPU 圖形處理器［10－13］是專門負(fù)責(zé)圖形渲染的核心處理器，內(nèi)部有上百個計算核心。CUDA 是Nvidia公司開發(fā)的一種SIMT 并行編程模型，用于開發(fā)高度并行的線程級程序。MPI［14－15］是一種基于消息傳遞并行編程模型的標(biāo)準(zhǔn)規(guī)范，從程序設(shè)計的角度理解是一個并行函數(shù)庫，其具體實(shí)現(xiàn)包括MPICH，OPEN MPI，LAM，IBM MPL 等多個版本，最常用和穩(wěn)定的是MPICH 和OPENMPI。MPI是目前超大規(guī)模并行計算應(yīng)用最廣、效率最高的并行編程模型，幾乎被所有并行環(huán)境（共享和分布式存儲并行機(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ìn)行高性能集群計算。

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

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

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

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

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

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

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

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

圖2 基于MPI和CUDA 的偏移計算流程

3 算法測試與效果分析

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

表1 測試環(huán)境

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

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

表2 測試計算時間 ms

表3 測試計算加速比

4 結(jié)束語

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

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

［1］馬昭軍，唐建明.疊前時間偏移在三維轉(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 平臺上的實(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］劉國峰，劉洪，王秀閩，等.Kirchhoff積分疊前時間偏移的兩種走時計算及并行算法［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í)現(xiàn)［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］劉國峰，劉欽，李博，等.油氣勘探地震資料處理GPU／CPU 協(xié)同 并行計算［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 高性能計算之CUDA［M］.北京：中國水利電力出版社，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］都志輝.高性能計算之并行編程技術(shù)－MPI并行程序設(shè) 計［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