海量地震數(shù)據(jù)疊前逆時偏移的多GPU聯(lián)合并行計算策略

孔祥寧，張慧宇，劉守偉，李晶晶

（1.中國石油化工股份有限公司石油物探技術(shù)研究院，江蘇南京 211103；2.同濟大學(xué)，上海 200092；3.中國石油大學(xué)（華東），山東青島 266580）

隨著高性能計算機及其應(yīng)用技術(shù)的高速發(fā)展，疊前逆時深度偏移方法已經(jīng)成為對復(fù)雜地質(zhì)體進行成像現(xiàn)實可行的手段，該方法在成像精度上相對于其它偏移方法的優(yōu)勢也被普遍認(rèn)可［1］。疊前逆時偏移方法直接求解雙程微分波動方程，避免了單程波偏移方法中上、下行波的分離，對波動方程的近似較少，克服了偏移傾角的限制，可以有效地處理縱、橫向速度存在劇烈變化的地球介質(zhì)物性特征［2－3］。該技術(shù)除了具有成像精度高，不受介質(zhì)縱、橫向速度變化和高陡傾角構(gòu)造的影響等優(yōu)點外，甚至還可以利用一些傳統(tǒng)處理中作為噪聲壓制的特殊波場（如多次波、回轉(zhuǎn)波及棱柱波等）進行成像。

疊前逆時深度偏移并不是一個新的理論方法，早在20世紀(jì)七八十年代就有學(xué)者提出了逆時偏移的基本思想和實現(xiàn)方法［4－6］。雖然RTM 成像精度高，且對速度的適應(yīng)性強，但是在實際生產(chǎn)中的應(yīng)用規(guī)模還非常小，這主要是由于RTM 對計算量和存儲量的需求巨大等瓶頸問題不易解決［7－9］。隨著大型計算集群的快速發(fā)展，計算需求和成本對RTM 的限制逐漸減小，RTM 的工程化應(yīng)用已成為可能［10］。近年在高性能計算方面出現(xiàn)的GPU／CPU 協(xié)同并行計算技術(shù)，迅速成為高性能密集計算發(fā)展的趨勢。此種異構(gòu)集群平臺非常適合于RTM 等大規(guī)模數(shù)據(jù)的處理［11］，與CPU 集群相比，具有投入少，后期運行維護成本低等優(yōu)點。國內(nèi)外在基于CPU／GPU 的地震疊前逆時偏移計算方面的研究工作不斷取得突破性的成果。Araya－Polo等［12］首先用Cell處理器實現(xiàn)了三維疊前深度逆時偏移；Micikevicius［13］給出了利用GPU 實現(xiàn)高階有限差分的算法；Darren 等［10］在利用GPU實現(xiàn)逆時偏移方面做了有益的嘗試；李博等［14］提出了地震疊前逆時偏移算法的CPU／GPU 實施對策；趙磊等［11］根據(jù)CPU／GPU 異構(gòu)平臺計算能力非常強的特點，以犧牲計算量換取存儲量和I／O 次數(shù)的減少，實現(xiàn)地震疊前逆時偏移計算。

簡單地說，GPU／CPU 協(xié)同并行計算就是將GPU 和CPU 兩種不同架構(gòu)的處理器結(jié)合在一起，組成硬件上的協(xié)同并行模式，同時在應(yīng)用程序編寫上實現(xiàn)GPU 和CPU 軟件協(xié)同配合的并行計算。我們主要介紹采用GPU 多卡聯(lián)合策略解決GPU顯存不足的瓶頸，在GPU／CPU 異構(gòu)集群平臺上實現(xiàn)RTM 的協(xié)同并行計算。

1 基本原理

各向同性介質(zhì)下，RTM 需要直接求解的三維聲波方程為

式中：u（x，y，z，t）為地表記錄的壓力波場；v（x，y，z）為縱橫向可變的介質(zhì)速度。

在求解方程（1）的過程中需要利用合適的成像條件提取成像值。本文中逆時偏移的實現(xiàn)采用了方程（2）所示成像條件，該成像條件考慮了震源波場的照明補償，較直接的相關(guān)成像條件具有較優(yōu)的保幅特性，可為后續(xù)的AVO 等屬性分析提供更真實的地震信息。

式中：tmax為最大記錄時間；Ss（x，y，z，t）為正向外推的震源波場；Rs（x，y，z，t）為反向外推的記錄波場；I（x，y，z）為點（x，y，z）的成像結(jié)果；σ為一小常量，以保證方程（2）的穩(wěn)定性。

逆時偏移的實現(xiàn)流程如圖1所示。

圖1 逆時偏移技術(shù)流程

逆時偏移方法是通過雙程波波動方程在時間上對地震資料進行反向外推并結(jié)合成像條件實現(xiàn)偏移的，對全波場進行逆時外推，避免了上、下行波的分離，且不受傾角的限制，并能夠?qū)θ我鈨A斜構(gòu)造甚至回轉(zhuǎn)波進行成像［14］。逆時偏移在算法實現(xiàn)上主要包括3部分關(guān)鍵技術(shù)，即震源波場的正向延拓、接收點波場的反向延拓和恰當(dāng)?shù)某上駰l件。由于震源波場沿時間是正向的，疊前炮記錄波場沿時間是逆向的，要將這兩個波場做零延遲互相關(guān)，就必須要保存其中一個波場，這就是逆時偏移面臨的存儲量問題。震源波場和檢波點波場都根據(jù)方程（1）進行雙程外推，無論是利用有限差分法求解還是利用偽譜法求解都需要巨大的計算量，這就是逆時偏移面臨的計算量問題。

目前，為了解決巨大的存儲需求問題，工業(yè)界普遍使用的是Check－pointing技術(shù)［15］，這種方法的核心思想是將偏移過程沿時間分成若干片段，每次偏移其中的一段時間，這樣只需存儲當(dāng)前時間片段中的波場即可，而付出的代價是增加0.5倍的計算量。這種方法較有效地減少了存儲需求量，但存儲量仍然較大，而且I／O 量非常巨大，極大地影響了計算效率。為此，我們采用震源波場重構(gòu)的存儲策略，即先將震源波場正向外推到Tmax，外推過程中只存儲邊界波場信息，然后在檢波點波場進行逆推時作為邊界條件再重構(gòu)震源波場，再應(yīng)用成像條件。這種方法同樣增加了0.5倍的計算量，但存儲量和I／O 量均得到了大幅度的降低。

對于計算量太大這一瓶頸問題，我們引入了GPU／CPU 聯(lián)合并行的計算策略。無論是波場的正向延拓還是反向延拓，我們都采用了顯式有限差分的計算方法，該方法計算密度非常高，并行性非常好，適用于利用GPU 進行并行計算。通過高性能的GPU／CPU 異構(gòu)集群平臺，充分利用GPU 在科學(xué)計算方面的強大優(yōu)勢，建立了實用的疊前逆時深度偏移成像技術(shù)流程。

2 RTM 并行計算策略

地震疊前逆時深度偏移是典型的大計算吞吐、大數(shù)據(jù)吞吐的地震數(shù)據(jù)處理任務(wù)，須采用并行計算的策略來實現(xiàn)其實際應(yīng)用；而RTM 中所固有的可分解性和線性疊加性質(zhì)，使其具有良好的并行能力，這為我們實現(xiàn)RTM 的并行計算提供了條件。GPU／CPU 協(xié)同并行計算就是將GPU 和CPU 兩種不同架構(gòu)的處理器結(jié)合在一起，組成硬件上的協(xié)同并行模式，同時在應(yīng)用程序編寫上實現(xiàn)GPU 和CPU 軟件協(xié)同配合的并行計算。CPU 主要負(fù)責(zé)GPU 的控制、數(shù)據(jù)的準(zhǔn)備、數(shù)據(jù)在節(jié)點間的發(fā)送與接收等，即進行并行控制；GPU 主要進行RTM 中最耗時的波場外推計算及進行并行計算。有限差分法波場外推是典型的單指令多數(shù)據(jù)的計算，非常適合用GPU 處理。

設(shè)計并行算法首先要選擇合適的并行粒度。并行算法根據(jù)計算任務(wù)的大小可分為粗粒度并行算法、細(xì)粒度并行算法和介于二者之間的中粒度并行算法。并行粒度的選擇在對物理問題并行性分析的基礎(chǔ)上，必須兼顧通訊開銷、計算與通訊的重疊程度、負(fù)載均衡及容錯處理等并行程序設(shè)計中的關(guān)鍵問題。

2.1 CPU 上的并行設(shè)計

RTM 計算以一個炮集數(shù)據(jù)為基本計算單元，每一炮數(shù)據(jù)的偏移成像都是一個相對獨立的作業(yè)，計算時相互之間不需要或很少需要進行數(shù)據(jù)交換，因而具有很強的并行性，很適合作為一個獨立的并行作業(yè)。這是典型的SPMD 算法，即在不同結(jié)點上用相同的程序?qū)Σ煌臄?shù)據(jù)體進行計算。不同節(jié)點上的計算任務(wù)不是同時開始或結(jié)束，各個節(jié)點任務(wù)的開始與結(jié)束和其它節(jié)點任務(wù)的開始與結(jié)束無關(guān)，并且不同節(jié)點上相同程序不同數(shù)據(jù)的計算或通信也不是同時進行的。每個節(jié)點上每一份作業(yè)做完，即得到一份新的作業(yè)，便又立即開始執(zhí)行，這大大增加了計算相對于通信的比重，有效地提高了計算效率?？紤]到計算節(jié)點間性能的差異及多用戶運行環(huán)境下最大限度地發(fā)揮各個節(jié)點的計算效率，我們采用了“任務(wù)池”分配方式（圖2）的主從并行計算模式來保證RTM 計算時的動態(tài)負(fù)載均衡。

具體的實現(xiàn)過程是：首先，主進程讀取計算參數(shù)，計算出RTM 成像的總炮數(shù)及相關(guān)參數(shù)，形成“任務(wù)池”；然后，各個從進程到主進程取得所要計算的“任務(wù)”，進行RTM 成像計算，將每一炮的成像結(jié)果保存在本地磁盤的臨時文件中，同時從主進程獲取新的“任務(wù)”，直到所有的“任務(wù)”完成；最后，收集各個從進程的計算結(jié)果，并輸出最終RTM 成像結(jié)果。

圖2 實現(xiàn)RTM 的主從并行計算模式

2.2 GPU 上的并行設(shè)計

GPU 上實現(xiàn)逆時偏移計算主要包括兩項關(guān)鍵技術(shù)，一是并行策略，二是存儲策略。

我們知道，GPU 上的并行計算屬于細(xì)粒度并行算法，它的并行結(jié)構(gòu)分為3個層次：線程、線程塊以及由線程塊組成的線程網(wǎng)格，它可以針對數(shù)據(jù)體中的每個元素進行并行計算，力度之細(xì)是CPU 無法比擬的。我們所實現(xiàn)逆時偏移的主要計算熱點為有限差分計算，有限差分算子可以抽象為向量乘法問題。因此，將整個RTM 計算過程中計算量最密集的波場延拓通過GPU 并行策略實現(xiàn)，最能提高計算效率。需要注意的是，利用高階有限差分法計算需要大量的內(nèi)存讀寫，以三維時間二階、空間8階差分網(wǎng)格為例，每計算一個網(wǎng)格點的值都需要讀取周圍25個網(wǎng)格點的數(shù)據(jù)，內(nèi)存讀取冗余度非常高。為此，我們考慮借用GPU 提供的片內(nèi)存儲器共享內(nèi)存來降低數(shù)據(jù)讀取延遲，提高計算效率［17］。目前最新的一款nVidia Tesla Kepler K20X 圖形處理器有2 688個核心，可以同時處理2 688個樣本，尤其是針對RTM 的單精度計算，具有非常高的處理效率。

我們主要采用的存儲策略是震源邊界波場存儲與重構(gòu)。在波場正傳時，記錄每一個計算時間步的震源波場的6個面，當(dāng)記錄波場逆向外推時再由波動方程及其邊值條件、初值條件重構(gòu)空間震源波場。其公式表示為

這樣節(jié)省了大量的存儲空間，降低了數(shù)據(jù)的通訊量，提高了計算效率。

此外，由于GPU 的顯存有限，當(dāng)計算一炮數(shù)據(jù)所需空間超過GPU 一塊卡的顯存容量時，則無法進行逆時偏移計算。為打破這一限制，我們采用了GPU 多卡聯(lián)合計算模式，即將GPU 多卡進行捆綁，共同進行單炮RTM 計算，如圖3所示。

為了使算法滿足任意大規(guī)模地震數(shù)據(jù)的處理，我們設(shè)計了利用內(nèi)存作為中轉(zhuǎn)站的數(shù)據(jù)交換機制，這種機制不限定每個節(jié)點GPU 卡的數(shù)量，根據(jù)計算需求動態(tài)規(guī)劃每炮計算需要的節(jié)點數(shù)和GPU卡的數(shù)量。這樣，一次數(shù)據(jù)交換，就需要4次GPU與CPU 之間的數(shù)據(jù)傳輸以及一次CPU 內(nèi)存間的傳輸，其中CPU 與GPU 的數(shù)據(jù)I／O 具有一定的訪存延遲。

圖3 GPU 多卡聯(lián)合計算模式

通過GPU 多卡處理可以解決GPU 顯存不足的限制，但另一方面引入了通訊量的問題。利用多卡聯(lián)合計算一炮數(shù)據(jù)，意味著將一個炮數(shù)據(jù)體分別用不同的GPU 卡計算，那么在計算過程中就需要卡與卡之間的數(shù)據(jù)交換，從而增加了GPU 卡之間的I／O 通訊量。為了優(yōu)化GPU 卡間的數(shù)據(jù)通訊，我們采用數(shù)據(jù)I／O 隱藏策略，將數(shù)據(jù)體按照GPU顯存進行劃分，每塊卡上的又可分為獨立計算部分和邊界的重疊計算部分。計算過程分為兩步進行（圖4）：①對重疊數(shù)據(jù)部分進行計算；②在計算獨立數(shù)據(jù)部分的同時，進行卡與卡之間重疊數(shù)據(jù)的交換。通過數(shù)據(jù)計算與數(shù)據(jù)I／O 通訊的同時進行，實現(xiàn)了數(shù)據(jù)I／O 的隱藏策略。在新一代的GPU架構(gòu)中（Fermi架構(gòu)）增加了GPUDirectTM技術(shù)，該技術(shù)可以直接讀取和寫入CUDA 主機內(nèi)存，消除不必要的系統(tǒng)內(nèi)存拷貝和CPU 開銷，還支持GPU之間以及類似NUMA 結(jié)構(gòu)的GPU 到其他GPU內(nèi)存直接訪問的P2P 的DMA 直接傳輸。這些功能為未來版本的GPU 和其它設(shè)備之間的直接點對點通信奠定了基礎(chǔ)。

圖4 數(shù)據(jù)I／O 的隱藏策略

3 應(yīng)用實例

我們應(yīng)用某海上三維地震數(shù)據(jù)對上述方法進行測試。該地震數(shù)據(jù)采集采用雙震源激發(fā)，共6條接收線，最大接收道數(shù)408，最大偏移距5 100m，道距12.5m，最大反射長度6 144ms，采樣間隔2ms。處理面元12.5m×12.5m。

本次試處理的目標(biāo)成像范圍約84km2，為保證足夠的偏移孔徑，資料輸入面積約為185km2，一共有45 061炮，數(shù)據(jù)量367G。偏移網(wǎng)格為：Ny＝997（dy＝12.5m），Nx＝1 188（dx＝12.5m），Nz＝601（dz＝8m）；偏移孔徑為xapert＝7 000m，yapert＝7 000m；采用子波為雷克子波，主頻30Hz；時間延拓長度為6.0s，步長0.5ms。在CPU 平臺進行疊前逆時深度偏移的單炮用時約為26.2h，而在GPU平臺實現(xiàn)疊前逆時深度偏移的單炮用時約為0.43h。圖5為基于CPU 平臺和GPU 平臺的逆時偏移剖面對比結(jié)果，可以看出，兩者成像效果相當(dāng)，均能對巖丘邊界進行較好的成像，清晰刻畫出巖體與周圍地層的接觸關(guān)系。

圖5 某海上三維地震數(shù)據(jù)RTM 偏移剖面

在GPU／CPU 異構(gòu)平臺下，RTM 單炮的計算用時基本不超過0.47h，且主要集中在0.3～0.4h（圖6）。試處理結(jié)果表明，利用CPU 實現(xiàn)逆時偏移耗時巨大，而利用GPU 計算則大大減少了計算周期，取得了60倍的加速比。因此，在GPU／CPU 平臺及多GPU 聯(lián)合并行計算策略的支持下，三維疊前逆時深度偏移已經(jīng)可以實現(xiàn)海量地震數(shù)據(jù)的實際生產(chǎn)性處理。

圖6 基于GPU／CPU 異構(gòu)平臺的RTM 每炮計算時間

4 結(jié)束語

通過震源波場重構(gòu)策略及GPU／CPU 協(xié)同并行計算，在保證成像品質(zhì)的前提下，充分發(fā)揮了GPU 的計算優(yōu)勢，極大地提高了RTM 的計算效率。而GPU 多卡聯(lián)合處理技術(shù)能夠有效地解決GPU 顯存不足這一瓶頸。實際資料的處理結(jié)果表明，借助于GPU／CPU 異構(gòu)集群后，RTM 這種高精度成像方法能夠達(dá)到實際應(yīng)用的要求，有利于對復(fù)雜地質(zhì)體進行準(zhǔn)確、有效地成像。

［1］劉定進，楊勤勇，方伍寶，等.疊前逆時深度偏移成像的實現(xiàn)與應(yīng)用［J］.石油物探，2011，50（6）：545－549 Liu D J，Yang Q Y，F(xiàn)ang W B，et al.Realization and practices of pre－stack reverse－time depth migration［J］.Geophysical Prospecting for Petroleum，2011，50（6）：545－549

［2］陳可洋.基于高階有限差分的波動方程疊前逆時偏移方法［J］.石油物探，2009，48（5）：475－478 Chen K Y.Wave equation pre－stack reverse－time migration method based on high－order finite－difference approach［J］.Geophysical Prospecting for Petroleum，2009，48（5）：475－478

［3］楊勤勇，段心標(biāo).逆時偏移技術(shù)發(fā)展現(xiàn)狀與趨勢［J］.石油物探，2010，49（1）：92－98 Yang Q Y，Duan X B.Development status and trend of reverse time migration［J］.Geophysical Prospecting for Petroleum，2010，49（1）：92－98

［4］Whitmore N D.Iterative depth migration by backward time propagation［J］.Expanded Abstracts of 53rdAnnual Internat SEG Mtg，1983，827－830

［5］McMechan G A.Migration by extrapolation of timedependent boundary values［J］.Geophysical Prospecting，1983，31（3）：413－420

［6］Loewenthal D，Mufti I R.Reversed time migration in spatial frequency domain［J］.Geophysics，1983，48（5）：627－635

［7］陳可洋.地震波逆時偏移方法研究綜述［J］.勘探地球物理進展，2010，33（3）：153－159 Chen K Y.Review of seismic reverse time migration methods［J］.Progress in Exploration Geophysics，2010，33（3）：153－159

［8］Yoon K，Marfurt K J，Starr W.Challenges in re－verse－time migration［J］.Expanded Abstracts of 74thAnnual Internat SEG Mtg，2004，1057－1060

［9］Dussaud E.Computational strategies for reversetime migration［J］.Expanded Abstracts of 78thAnnual Internat SEG Mtg，2008，2267－2271

［10］Foltinek D，Eaton D，Mahovsky J，et a1.Industrialscale reverse time migration on GPU hardware［J］.Expanded Abstracts of 79thAnnual Internat SEG Mtg，2009，2789－2793

［11］趙磊，王華忠，劉守偉.逆時深度偏移成像方法及其在CPU／GPU 異構(gòu)平臺上的實現(xiàn)［J］.巖性油氣藏，2010，22（F07）：36－41 Zhao L，Wang H Z，Liu S W.Reverse－time depth migration imaging and its application at CPU／GPU platform［J］.Lithologic Reservoirs，2010，22（F07）：36－41

［12］Araya－Polo M，Rubio F，de la Cruz R，et a1.3Dseismic imaging through reverse－time migration on homogeneous and heterogeneous multi－core processors［J］.Journal of Scientific Programming，2009，17（1／2）：186－198

［13］Micikevicius P.3Dfinite difference computation on GPUs using CUDA［J］.Proceedings of 2ndWorkshop on General Purpose Processing on Graphics Processing Units，2009，79－84

［14］李博，劉紅偉，劉國峰，等.地震疊前逆時偏移算法的CPU／GPU實施對策［J］.地球物理學(xué)報，2010，53（12）：2938－2943 Li B，Liu H W，Liu G F，et al.Computational strategy of seismic pre－stack reverse time migration on CPU／GPU［J］.Chinese Journal of Geophysics，2010，53（12）：2938－2943

［15］Symes W W.Reverse time migration with optimal check pointing［J］.Geophysics，2007，72（5）：SM213－SM221

［16］Villarreal A，Scales J A.Distributed three dimensional finite－difference modeling of wave propagation in acoustic media［J］.Computers in Physics，1997，11（4）：388－399

［17］李博，劉國峰，劉洪.地震疊前時間偏移的一種圖形處理器提速實現(xiàn)方法［J］.地球物理學(xué)報，2009，52（1）：245－252 Li B，Liu G F，Liu H.A method of using GPU to accelerate seismic pre－stack time migration［J］.Chinese Journal of Geophysics，2009，52（1）：245－252