沈驥千

（江蘇煤炭地質(zhì)物測隊，南京 210046）

三維彈性波逆時偏移中多卡GPU的應(yīng)用

沈驥千

（江蘇煤炭地質(zhì)物測隊，南京 210046）

地震勘探對象的日益復(fù)雜化，對成像效果提出了更高的要求。逆時偏移在理論上可以對復(fù)雜構(gòu)造準(zhǔn)確成像，且二維逆時偏移的研究也取得了一系列成果。而逆時偏移巨大的存儲、計算量以及野外單炮數(shù)據(jù)量的龐大，已成為制約其在三維多分量資料中應(yīng)用的瓶頸。針對有效邊界存儲在三維逆時偏移中存儲較大的問題，采用隨機(jī)邊界條件以降低逆時偏移存儲量并提高偏移效率。該邊界條件無需保存各分量在邊界處的波場值，即可實現(xiàn)正向延拓波場的逆時重構(gòu)；考慮到三維逆時偏移需開辟較大數(shù)組空間用于計算，提出應(yīng)用區(qū)域分解技術(shù)將計算數(shù)據(jù)分配到不同節(jié)點(diǎn)，實現(xiàn)基于MPI+CUDA的協(xié)同并行。SEG/EAGE鹽丘模型算例證明了上述方法的有效性。

三維彈性波逆時偏移；隨機(jī)邊界；區(qū)域分解；MPI+CUDA協(xié)同并行

0 引言

基于標(biāo)量波理論的反射縱波勘探技術(shù)在以往的能源開發(fā)過程中發(fā)揮了重要作用，但隨著油氣開發(fā)對地震勘探精度要求的不斷提高和勘探目標(biāo)構(gòu)造和巖性復(fù)雜程度的不斷增加，縱波勘探精度的提高越來越受制于其理論上的局限與波場信息單一等缺陷，在碳酸鹽巖裂縫性儲層、煤層氣和頁巖氣儲層等的勘探方面難以獲得滿意效果。近年來，以矢量波動理論為基礎(chǔ)的多波地震勘探技術(shù)引起了業(yè)界的廣泛重視，同常規(guī)縱波勘探技術(shù)相比，多波地震技術(shù)的理論假設(shè)與實際情況更為接近，且橫波在復(fù)雜介質(zhì)中具有獨(dú)特的傳播機(jī)理，故理論上多波地震技術(shù)解決實際問題的能力更強(qiáng)，更有利于提高勘探精度或?qū)崿F(xiàn)某些特殊目標(biāo)的精確勘探。

疊前深度偏移是多波地震技術(shù)的重要研究內(nèi)容，對多波地震資料進(jìn)行疊前深度偏移處理的意義在于：一方面，地下構(gòu)造的精確成像需要借助高精度的疊前偏移技術(shù)來實現(xiàn)；另一方面，多波偏移剖面中所包含的地震波動力學(xué)特征還可用于巖性和流體識別；同時偏移得到的共成像點(diǎn)道集還可為振幅隨偏移距變化（AVO）、振幅隨入射角變化（AVA）或其它疊前反演工作提供輸入數(shù)據(jù)。因此，多波地震資料偏移應(yīng)當(dāng)實現(xiàn)兩個目標(biāo)：（1）地下構(gòu)造的縱、橫波準(zhǔn)確成像；（2）獲取保真的縱橫波疊前反演道集和疊加剖面。

基于彈性波理論的多分量地震資料聯(lián)合逆時偏移技術(shù)是實現(xiàn)多波地震縱橫波偏移成像的有效工具[1]。一方面，彈性波逆時偏移技術(shù)以雙程波動方程為理論出發(fā)點(diǎn)，它無傾角限制，理論上可實現(xiàn)任意復(fù)雜構(gòu)造的成像[2-5]；另一方面，相較于標(biāo)量波動方程逆時偏移，多分量聯(lián)合逆時偏移保留了地震波的彈性、矢量特征，可以對氣云帶等特殊目標(biāo)獲得良好成像效果[6]。

逆時偏移思想提出之初，在疊后地震資料偏移成像中獲得了良好的成像剖面[3]；但受限于計算機(jī)設(shè)備的發(fā)展，其巨大的計算量和存儲I/O（輸入/輸出）成為其工業(yè)化應(yīng)用的巨大瓶頸。隨著計算機(jī)硬件設(shè)備的革新和GPU（Graphic Processing Unit）以及其編程架構(gòu)CUDA（Computing Unified Device Archi?tecture）的快速發(fā)展，標(biāo)量波動方程逆時偏移并行化算法已漸趨成熟[7-10]。但多分量地震資料聯(lián)合逆時偏移的并行加速尚存在諸多問題，主要表現(xiàn)在：①多分量地震聯(lián)合逆時偏移存儲量和計算量遠(yuǎn)大于標(biāo)量波動方程逆時偏移；②基于CPU/GPU協(xié)同并行的標(biāo)量波動方程逆時偏移相關(guān)配套技術(shù)應(yīng)用于彈性波逆時偏移尚需改進(jìn)，如有效邊界存儲等。

目前的野外采集普遍向高密度、三維方向發(fā)展，為確保復(fù)雜構(gòu)造觀測波場的連續(xù)性，多采用小炮距、小接收線距、長排列接收，部分探區(qū)單炮記錄過萬道[11]。這導(dǎo)致三維彈性波逆時偏移過程中開辟的數(shù)組空間所占內(nèi)存達(dá)10 G，依靠單節(jié)點(diǎn)GPU的加速不能滿足存儲需求。另一方面，在二維波動方程逆時偏移中為解決逆時延拓與正向延拓時間序列不一致性，常采用有效邊界[12]存儲正向延拓邊界波場，實現(xiàn)正演波場的重構(gòu)來保證與逆時延拓的時間一致性。但是在三維多分量地震資料聯(lián)合逆時偏移中，需要保存的邊界波場比之于二維情況增加數(shù)個數(shù)量級，加重了存儲負(fù)擔(dān)和存取I/O。

本文針對有效邊界存儲在三維情況下存儲過大的情況，提出應(yīng)用隨機(jī)邊界條件，在不保存邊界波場的情況下，實現(xiàn)正向延拓波場的重構(gòu)；考慮三維逆時偏移情況下，開辟數(shù)組空間內(nèi)存較大，提出應(yīng)用區(qū)域分解技術(shù)將計算數(shù)據(jù)分配到不同節(jié)點(diǎn)，實現(xiàn)基于MPI+CUDA的協(xié)同并行。

1 彈性波逆時偏移的基本原理

彈性波逆時偏移將多分量地震記錄當(dāng)作彈性波場的邊值問題[13]來處理，由計算機(jī)實現(xiàn)炮點(diǎn)波場的正向延拓和接收點(diǎn)波場的逆時延拓，并在延拓過程中采用相應(yīng)的算法實現(xiàn)縱橫波的成像。其輸入數(shù)據(jù)一般為未經(jīng)坐標(biāo)旋轉(zhuǎn)的三分量炮集，輸出數(shù)據(jù)一般包括反射縱波成像結(jié)果和轉(zhuǎn)換橫波成像結(jié)果。

彈性波逆時偏移的基本原理為：假定地下介質(zhì)由一系列繞射點(diǎn)組成，地面、海底或井中接收到的三分量記錄為各繞射點(diǎn)產(chǎn)生的縱橫繞射波場在接收點(diǎn)處的疊加響應(yīng)，其中各繞射點(diǎn)縱橫繞射波的產(chǎn)生時間為炮點(diǎn)波場的主能量到達(dá)該點(diǎn)的時間[14]，繞射點(diǎn)產(chǎn)生的縱橫繞射波振幅取決于該點(diǎn)的反射系數(shù)[15]，由此可通過求取各繞射點(diǎn)產(chǎn)生縱橫繞射波的時間和該時刻的縱橫波振幅來確定地下各繞射點(diǎn)的空間位置和反射系數(shù)，實現(xiàn)地下介質(zhì)的縱橫波成像。實際成像過程一般由以下幾部分組成：①炮點(diǎn)波場的正向延拓和重構(gòu)；②接收點(diǎn)波場的逆時延拓；③炮、檢波場延拓過程中的縱橫波解耦與橫波的標(biāo)量化處理；④縱橫波成像。

1.1 一階速度——應(yīng)力彈性波方程及其高階延拓差分格式

根據(jù)各向同性介質(zhì)中，位移-應(yīng)變關(guān)系、應(yīng)力-應(yīng)變關(guān)系、應(yīng)力-位移關(guān)系[16]可推導(dǎo)得到各向同性介質(zhì)中的一階應(yīng)力-速度彈性波方程：

其中，υx、υy、υz為質(zhì)點(diǎn)振動速度；σxx、σyy、σzz、σxy、σxz、σyz為應(yīng)力分量；ρ為密度；υp為縱波速度，υs為橫波速度；x、y、z分別為空間坐標(biāo)；t為時間。

在交錯網(wǎng)格空間中對式（1）進(jìn)行差分離散可得到時間2階，空間任意偶數(shù)階的差分格式，以式（1）中第一個方程為例：

其中，Δx、Δy、Δz分別表示坐標(biāo)系中x、y、z方向上的網(wǎng)格間距；Δt表示時間采樣間隔；i，j，k分別代表x、y、z方向上離散點(diǎn)的序號；n表示時間上的離散點(diǎn)的序號。且ρ、μ半網(wǎng)格點(diǎn)上的值由同一方向上相鄰整網(wǎng)格點(diǎn)求取算術(shù)平均得到。2N為差分階數(shù)，D(mN)為2N階交錯網(wǎng)格差分系數(shù)，其計算公式為[16]：

1.2 穩(wěn)定性條件與吸收邊界條件

彈性波方程交錯網(wǎng)格有限差分格式屬于顯示差分型遞推格式，計算當(dāng)前時刻波場值時，需利用上一時刻波場值。由于數(shù)值計算截斷誤差的影響，隨著時間遞推誤差會累積，破壞數(shù)值穩(wěn)定性甚至引起計算溢出[17]。

通過對三維各向同性介質(zhì)一階雙曲型彈性波方程進(jìn)行平面諧波分析[16]可得穩(wěn)定性條件為：

式（4）中，Δx、Δy、Δz分別表示坐標(biāo)系中x、y、z方向上的網(wǎng)格間距；Δt表示時間采樣間隔；υmax為背景縱、橫波速度的極大值。對于三維各向同性彈性介質(zhì)，地震波傳播速度與傳播方向無關(guān)，所以，在三個方向上的穩(wěn)定條件相同。

采用無分裂式PML邊界條件解決階段邊界的偽反射問題[18]，依據(jù)無分裂PML邊界條件的實現(xiàn)思路，可得到式（1）基于無分裂PML邊界條件的控制方程；以式（1）中第一個方程為例：

上式中Ωxx、Ωxy、Ωxz為引入的中間變量，其所對應(yīng)的控制方程為：

其中，d（x）、d（y）、d（z）為三個方向的衰減因子；根據(jù)公式（5）、（6）可依此類推基于無分裂PML邊界條件的彈性波控制方程，其計算步驟為：（1）對包含PML邊界的整個區(qū)域進(jìn)行常規(guī)的彈性波波場求解；（2）在PML邊界區(qū)域?qū)⒌玫降某Ｒ?guī)波場減去衰減量，然后按照時間軸方向更新輔助變量。由控制方程和計算步驟可知，無分裂的PML邊界條件避免了對彈性波波場進(jìn)行分裂，降低了存儲。由于采用同一套方程也不會出現(xiàn)在邊界處兩套方程耦合的編程難題。

1.3 互相關(guān)成像條件

成像條件決定著逆時偏移成像剖面的質(zhì)量。目前業(yè)界廣泛應(yīng)用互相關(guān)成像條件，其基本原理為：利用炮點(diǎn)波場和接收點(diǎn)逆時延拓波場做零延時互相關(guān)，在反射點(diǎn)位置的波場將得到加強(qiáng)，非反射點(diǎn)處的波場將削弱，最后所有時刻波場相加得到偏移剖面。互相關(guān)成像條件的實現(xiàn)算法為：

其中，Ipp（x，z）為PP波偏移剖面，Ips（x，z）為PS波偏移剖面，Isp（x，z）為SP波偏移剖面，Iss（x，z）為SS波偏移剖面，S_P（x，z）、S_S（x，z）為炮點(diǎn)縱、橫波波場，R_P（x，z）、R_S（x，z）為接收點(diǎn)縱波和橫波波場，純縱橫波波場的獲得可參考文獻(xiàn)[19，20]，tmax為地震記錄長度。

2 基于多卡GPU的三維彈性波逆時偏移及其相關(guān)技術(shù)

上述彈性波逆時偏移基本流程中，在利用互相關(guān)成像條件時需同時應(yīng)用相同時刻的檢波點(diǎn)波場和炮點(diǎn)波場，上述兩個波場的時間序列相反，因此需要重構(gòu)炮點(diǎn)波場以實現(xiàn)縱橫波成像；本文采用隨機(jī)邊界條件散射人工邊界的偽反射，實現(xiàn)波場的重構(gòu)。為了解決三維情況下單卡GPU顯存不足的難題，利用區(qū)域分解技術(shù)將數(shù)據(jù)劃分，實現(xiàn)大規(guī)模數(shù)據(jù)的三維彈性波逆時偏移。

2.1 隨機(jī)邊界條件的應(yīng)用

在解決炮點(diǎn)、檢波點(diǎn)波場時間序列不一致的問題中，國內(nèi)外學(xué)者先后提出Checkpointing方法、有效邊界存儲策略[12]，來降低存儲量。在二維彈性波逆時偏移中，有效邊界存儲策略已獲得業(yè)界廣泛認(rèn)可；有效邊界存儲只需保存PML區(qū)域差分階數(shù)層的邊界值，即可實現(xiàn)波場的準(zhǔn)確重構(gòu)。但在三維情況下，其存存儲量仍相當(dāng)可觀，假設(shè)三維模型網(wǎng)格點(diǎn)個數(shù)為200 m×300 m×400 m，采樣點(diǎn)數(shù)為5 000，有限差分階數(shù)為12階，若采用有效邊界存儲其存儲量達(dá)60G左右，這顯然難以滿足生成要求。因此，提出使用隨機(jī)邊界對波場進(jìn)行散射[21]，因隨機(jī)邊界不涉及波場的衰減，其過程可逆，學(xué)者們提出“以計算換存儲”，利用最后時刻波場即可重構(gòu)得到每一時刻震源波場值。

隨機(jī)邊界條件實質(zhì)為在計算區(qū)域外部隨機(jī)賦速度值，利用隨機(jī)速度的無規(guī)律性，將波場值散射，這樣在進(jìn)行互相關(guān)成像和疊加時，無規(guī)律的波場值能量得到減弱，反射點(diǎn)成像值相對得到增強(qiáng)。所以，隨機(jī)邊界條件的關(guān)鍵之處在于隨機(jī)邊界層內(nèi)隨機(jī)速度的選擇，常用隨機(jī)速度賦值公式為：

式中，Vrandom（x，y，z）為邊界層內(nèi)（x，y，z）處應(yīng)取的隨機(jī)速度，v（x，y，z）為對應(yīng)的背景速度，rand是隨機(jī)函數(shù)求得的隨機(jī)數(shù)，d為隨機(jī)邊界層內(nèi)的點(diǎn)（x，y，z）與區(qū)域邊界的距離，v0（x，y，z）為最大穩(wěn)定速度。

為了驗證隨機(jī)邊界條件的適用性，本文采用二維均勻速度模型，驗證隨機(jī)邊界的重構(gòu)和散射效果。二維均勻速度模型大小1500 m×1500 m，縱波速度為2500 m/s，橫波速度1800 m/s，炮點(diǎn)位置位于（750 m，750 m），縱波源激發(fā)，記錄不同時刻vz分量的正演和重構(gòu)波場快照。圖1即為均勻介質(zhì)隨機(jī)邊界速度場，圖2、圖3、圖4即為不同時刻正向延拓和逆時重構(gòu)的波場值，由此可驗證隨機(jī)邊界的有效性及波場重構(gòu)的正確性。

2.2 區(qū)域分解技術(shù)及基于MPI+CUDA的彈性波逆時偏移流程

圖1 二維均勻介質(zhì)隨機(jī)邊界速度場Figure1 2D homogeneous medium stochastic boundary velocity field

圖2 二維均勻介質(zhì)隨機(jī)邊界正向延拓（250 ms）及重構(gòu)波場快照Figure 2 2D homogeneous medium stochastic boundary forward continuation（250 ms）and wave field reconstruction snapshot

GPU中有大量晶體管，是專為執(zhí)行復(fù)雜的數(shù)學(xué)和幾何計算而設(shè)計的。當(dāng)今的GPU已不再將這些計算局限于圖形渲染，其通用計算技術(shù)的發(fā)展和天生的并行處理機(jī)制在計算方面提供了高于CPU數(shù)十倍乃至上百倍于CPU的性能。業(yè)界已經(jīng)利用GPU的并行特性在二維逆時偏移取得了一系列突破[22-24]，然而對于大規(guī)模數(shù)據(jù)處理通常需要較大的計算機(jī)內(nèi)存（顯存），GPU目前的顯存一般為6～8 G，對于模型較大和實際數(shù)據(jù)的三維逆時偏移這顯然是不能滿足計算要求，需要進(jìn)行特殊情況考慮。本文的解決策略為將地震數(shù)據(jù)進(jìn)行區(qū)域分解[10]，形成多GPU協(xié)同處理單炮地震數(shù)據(jù)的快速計算方法；此方法在理論上可以對任意規(guī)模的三維地震資料進(jìn)行逆時偏移成像。

圖3 二維均勻介質(zhì)隨機(jī)邊界正向延拓（300 ms）及重構(gòu)波場快照Figure 3 2D homogeneous medium stochastic boundary forward continuation（300 ms）and wave field reconstruction snapshot

圖4 二維均勻介質(zhì)隨機(jī)邊界正向延拓（350 ms）及重構(gòu)波場快照Figure 4 2D homogeneous medium stochastic boundary forward continuation（350 ms）and wave field reconstruction snapshot

在對彈性波方程的差分格式，進(jìn)行有限差分計算編程時，將一個維度上（以y方向為例）的區(qū)域進(jìn)行劃分，按照實際資料的大小和顯卡顯存的限制來決定劃分的數(shù)量。利用區(qū)域分解進(jìn)行差分計算，需要解決的核心問題在于如何隱藏卡與卡間的數(shù)據(jù)交換。如圖5所示，為了保證差分計算的正確性，在每一個時間步需要在虛線區(qū)域進(jìn)行數(shù)據(jù)交換。

圖5 三維數(shù)據(jù)體區(qū)域分解示意圖Figure 5 A schematic diagram of 3D data volume domain decomposition

以其中的兩塊卡為例，說明分解區(qū)域后每塊卡的計算與數(shù)據(jù)傳遞。如圖6所示，每塊卡可以分為三部分（圖中分別用A、B和C表示），C區(qū)域為兩塊卡之間的重疊部分，厚度由差分格式來決定，B區(qū)域為傳遞部分與C區(qū)域保持一樣的厚度，A區(qū)域為顯卡本身的計算區(qū)域。在正演和延拓計算的每一個時間步內(nèi)，C區(qū)域的網(wǎng)格處在邊界而無法進(jìn)行差分運(yùn)算，所以需要將B區(qū)域的數(shù)據(jù)傳遞到C。這樣計算過程就分為了如下步驟：首先是每塊卡計算B區(qū)域；然后，每塊卡計算A區(qū)域的同時，傳遞B區(qū)域數(shù)據(jù)到相鄰顯卡。

圖6 區(qū)域分解數(shù)據(jù)交換示意圖Figure 6 A schematic diagram of domain decomposition data exchange

計算能力高于2.0的NVIDIA顯卡支持在同節(jié)點(diǎn)上進(jìn)行顯卡與顯卡之間的數(shù)據(jù)交換。而本文提出的多卡并行是基于多節(jié)點(diǎn)的并行，在數(shù)據(jù)傳遞時需要將傳遞數(shù)據(jù)通過內(nèi)存中轉(zhuǎn)傳遞，由MPI的每一個進(jìn)程控制每塊顯卡，在MPI進(jìn)程組內(nèi)進(jìn)行非阻塞的數(shù)據(jù)交換。整個的實現(xiàn)流程圖如圖7所示：①由主進(jìn)程讀取到炮記錄和速度模型并分發(fā)于進(jìn)程；②各進(jìn)程控制GPU顯卡實現(xiàn)每一時間步的圖4B區(qū)域波場計算；③GPU回傳數(shù)據(jù)并在內(nèi)存中交換，同時實現(xiàn)圖4A區(qū)域的波場計算；④循環(huán)迭代波場，最后由主進(jìn)程規(guī)約得到最終的偏移結(jié)果。

3 數(shù)值算例

SEG/EAEG鹽丘模型是一個國際上通用測試三維復(fù)雜構(gòu)造成像效果的地質(zhì)模型。本文采用SEG/ EAEG鹽丘速度模型作為縱波模型，通過固定的縱橫波速度比構(gòu)造橫波模型，炮集記錄由一階彈性波方程有限差分正演獲得。模型參數(shù)如下：模型大小3380 m×3380 m×1050 m，網(wǎng)格大小Δx=Δy=Δz=5 m，時間采樣間隔0.35 ms，記錄長度2.8 s，震源為35 Hz雷克子波。觀測系統(tǒng)參數(shù)為：總炮數(shù)4690炮，67條測線，炮線間距50 m，炮間距60 m，接收線間距10 m，接收道間距5 m，每條測線70炮，每炮30400道接收；偏移算法采用本文所提出的基于MPI+CUDA的三維彈性波逆時偏移并行算法。圖8為三維鹽丘模型的縱橫波速度模型以及多分量逆時偏移剖面。

另外，PP、PS的成像效果較SP、SS的成像效果更好。原因為：PP，PS波是在縱波傾斜入射到彈性界面時產(chǎn)生能量較強(qiáng)，而SP、SS是由縱波傾斜入射到彈性界面產(chǎn)生的轉(zhuǎn)換橫波再次入射得到的地震記錄，該情況涉及波型的多次轉(zhuǎn)換，情況復(fù)雜。因此PP成像結(jié)果主要代表了P波的反射系數(shù)，PS主要成像結(jié)果代表了PS波的反射系數(shù)；而SP、SS成像結(jié)果并不對應(yīng)一種波的反射系數(shù)，多次轉(zhuǎn)換后并不清楚它的能量分布。PP波和PS波均可以對鹽丘模型構(gòu)造得到準(zhǔn)確的偏移成像結(jié)果。PP波主頻比PS波主頻低，其原因在于：模型資料的縱橫波記錄是由同一縱波源激發(fā)得到，且轉(zhuǎn)換橫波在空間域中的波長小于縱波波長，所以在模型資料的深度偏移成像結(jié)果中PS成像結(jié)果分辨率高于PP成像結(jié)果。

圖7 三維彈性波逆時偏移協(xié)同并行流程圖Figure 7 3D elastic wave RTM cooperative parallel flow chart

4 結(jié)論

隨機(jī)邊界條件在不保存邊界波場的情況下，可以實現(xiàn)正向延拓波場的逆時重構(gòu)，避免了三維彈性波逆時偏移大規(guī)模的邊界存儲問題；區(qū)域分解技術(shù)的應(yīng)用使得大規(guī)模數(shù)據(jù)資料的逆時偏移處理成為可能；在上述基礎(chǔ)上，本文提出的基于MPI+CUDA協(xié)同并行的三維彈性波逆時偏移方法可以對模型資料準(zhǔn)確成像。

[1]Gaiser J，Strudley A.Acquisition and application of multicomponent vector wavefields：are they practical[J].first break，2005，23(6)：61-67.

[2]Whitmore，N.D.Iterative depth migration by backward time propaga?tion.1983 SEG Annual Meeting.Society of Exploration Geophysicists，1983.

[3]McMechan G A.Migration by extrapolation of time-dependent bound?ary values.Geophysical Prospecting，1983，31(3):413-420.

[4]Baysal E，Kosloff D D，Sherwood J W C.Reverse time migration.Geo?physics，1983，48(11):1514-1524.

[5]Loewenthal D，Mufti I R.Reverse time migration in the spatial fre?quency domain.Geophysics，1983，48(5):627-635.

[6]楊佳佳.多分量地震波逆時偏移的關(guān)鍵技術(shù)研究[D].中國海洋大學(xué)，2015.

[7]何兵壽，張會星，韓月.雙程聲波方程疊前逆時深度偏移及其并行算法.煤炭學(xué)報，2010，35(3)：458-462.

[8]李博，劉國峰，劉洪.地震疊前逆時偏移的一種圖形處理器提速實現(xiàn)方法[J].地球物理學(xué)報，2009，52(12)：245-252.

[9]劉紅偉，李博，劉洪等.地震疊前逆時偏移高階有限差分算法及GPU實現(xiàn)[J].地球物理學(xué)報，2010，53(7)：1725-1733.

[10]劉守偉，王華忠，陳生昌，等.三維逆時偏移CPU/GPU機(jī)群實現(xiàn)方案研究[J].地球物理學(xué)報，2013，56(10).

[11]唐祥功，匡斌，杜繼修.多GPU協(xié)同三維疊前逆時偏移方法研究與應(yīng)用：石油地球物理勘探，2013，48（6）：910-914

[12]王保利，高靜懷，陳文超等.地震疊前逆時偏移的有效邊界存儲策略[J].地球物理學(xué)報，2012，55(7):2412-2421.

[13]Chang，W.F，McMechan G A.Elastic reverse-time migration[J].Geo?physics，1987，52(10)：1365-1375.

圖8 三維彈性波速度模型及逆時偏移成像結(jié)果Figure 8 3D elastic wave velocity model and RTM imaging result

[14]Claerbout，J.F.Toward a unified theory of reflector mapping[J].Geo?physics，1971，36(3):467-481.

[15]Zhi L，Chen S，Li X.Joint AVO inversion of PP and PS waves using exact Zoeppritz equation[M]//SEG Technical Program Expanded Ab?stracts 2013.Society of Exploration Geophysicists，2013：457-461.

[16]牟永光，裴正林.三維復(fù)雜介質(zhì)地震數(shù)值模擬[M].北京:石油工業(yè)出版社，2005，173-176.

[17]董良國，馬在田，曹景忠.一階彈性波方程交錯網(wǎng)格高階差分解法穩(wěn)定性研究[J].地球物理學(xué)報，2000，06:856-864.

[18]王鵬飛，何兵壽.聲波方程逆時偏移中的無分裂PML吸收邊界條件[J].工程地球物理學(xué)報，2015，12(5):583-590.

[19]Dellinger J，Etgen J.Wave-field separation in two-dimensional anisotropic media.Geophysics，1990，55(7)：914-919.

[20]Du Q Z，Gong X F，Zhang M，et al.3D PS-wave imaging with elastic reverse-time migration.Geophysics，2014，79(5)：S173-S184.

[21]Robert G.Clapp.Reverse time migration with random boundar?ies.79th Annual International Meeting，SEG Expanded Abstracts，2009，28：2809-2813.

[22]Micikevicius P，3D finite difference computation on GPUs using CU?DA.//Processing of 2th Workshop on General Purpose Processing on Graphics Processing Units，Expanded Abstracts，2008:79-84.

[23]Michea D，Komatitsch D，Mahovsky J，et al.Industrial-scale reverse time migration on GPU hardware.SEG Expand Abstracts，2009:2789-2793.

[24]Sun X Y，Suh S.Maxing throughput for high performance TTI-RTM: From CPU-GPU.SEG Expanded Abstracts，2011，182(1):389-402.

Application of Multi-card GPU in 3D Elastic Wave Reverse-time Migration(RTM)

Shen Jiqian
(Jiangsu Provincial Coal Geological Exploration,Geophysical Prospecting and Surveying Team,Nanjing,Jiangsu 210046)

With the increasing complexity of seismic exploration objects,higher requirements are put forward on imaging results.The RTM can be used to describe complex structures in theory.A series of research findings have been acquired in 2D RTM;however,the huge amount of storage,computational workload and field single-shot data in RTM have become the bottleneck to restrict its applica?tion in 3D multi-component data.In allusion to the issue of larger efficient boundary storage in 3D RTM,the stochastic boundary ap?proach is used to lower down RTM storage and improve migration efficiency.This boundary condition needs not to keep wave field val?ue of components at boundary and realize reverse-time reconstruction.Whereas 3D RTM needs larger array space in computation,the author applies the domain decomposition technique to distribute computational data to different nodes and realize cooperative parallel based on MPI+CUDA.The SEG/EAGE salt dome example has proved the effectiveness of the method.

3D elastic wave RTM;stochastic boundary;domain decomposition;MPI+CUDA cooperative parallel

P631.4

A

10.3969/j.issn.1674-1803.2017.03.14

1674-1803(2017)03-0065-07

國家自然科學(xué)基金（41674118），國家重大科技專項（2016ZX05027-002）聯(lián)合資助

沈驥千（1970—），男，高級工程師，主要從事地震勘探資料處理工作。

2017-02-10

責(zé)任編輯：孫常長