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

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

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

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

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

1 基本原理

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

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

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

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

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

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

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

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

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

2 RTM 并行計(jì)算策略

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 應(yīng)用實(shí)例

我們應(yīng)用某海上三維地震數(shù)據(jù)對(duì)上述方法進(jìn)行測(cè)試。該地震數(shù)據(jù)采集采用雙震源激發(fā)，共6條接收線，最大接收道數(shù)408，最大偏移距5 100m，道距12.5m，最大反射長(zhǎng)度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；時(shí)間延拓長(zhǎng)度為6.0s，步長(zhǎng)0.5ms。在CPU 平臺(tái)進(jìn)行疊前逆時(shí)深度偏移的單炮用時(shí)約為26.2h，而在GPU平臺(tái)實(shí)現(xiàn)疊前逆時(shí)深度偏移的單炮用時(shí)約為0.43h。圖5為基于CPU 平臺(tái)和GPU 平臺(tái)的逆時(shí)偏移剖面對(duì)比結(jié)果，可以看出，兩者成像效果相當(dāng)，均能對(duì)巖丘邊界進(jìn)行較好的成像，清晰刻畫出巖體與周圍地層的接觸關(guān)系。

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

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

圖6 基于GPU／CPU 異構(gòu)平臺(tái)的RTM 每炮計(jì)算時(shí)間

4 結(jié)束語(yǔ)

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

［1］劉定進(jìn)，楊勤勇，方伍寶，等.疊前逆時(shí)深度偏移成像的實(shí)現(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］陳可洋.基于高階有限差分的波動(dòng)方程疊前逆時(shí)偏移方法［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í)偏移技術(shù)發(fā)展現(xiàn)狀與趨勢(shì)［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］陳可洋.地震波逆時(shí)偏移方法研究綜述［J］.勘探地球物理進(jìn)展，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］趙磊，王華忠，劉守偉.逆時(shí)深度偏移成像方法及其在CPU／GPU 異構(gòu)平臺(tái)上的實(shí)現(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］李博，劉紅偉，劉國(guó)峰，等.地震疊前逆時(shí)偏移算法的CPU／GPU實(shí)施對(duì)策［J］.地球物理學(xué)報(bào)，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］李博，劉國(guó)峰，劉洪.地震疊前時(shí)間偏移的一種圖形處理器提速實(shí)現(xiàn)方法［J］.地球物理學(xué)報(bào)，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