碳酸鹽巖儲層疊前逆時偏移及GPU/CPU異構(gòu)計算方案研究

王闊，陳紅靈，郭聰，張文武，葉飛

（東北石油大學(xué)， 黑龍江 大慶 163318）

碳酸鹽巖儲層疊前逆時偏移及GPU/CPU異構(gòu)計算方案研究

王闊，陳紅靈，郭聰，張文武，葉飛

（東北石油大學(xué)， 黑龍江 大慶 163318）

針對碳酸鹽巖儲層地質(zhì)構(gòu)造的特點，提出采用基于聲波方程的疊前逆時偏移算法對其進行成像，并在偏移過程中引入優(yōu)化差分系數(shù)和完全匹配層(PML)邊界條件來提高波場的數(shù)值模擬精度，同時利用GPU/CPU異構(gòu)并行加速技術(shù)對算法進行優(yōu)化，從而顯著提升偏移的計算效率。模型測試表明，與傳統(tǒng)差分系數(shù)相比，同階的優(yōu)化差分系數(shù)具有更高精度，壓制頻散效果明顯；PML邊界條件能夠有效的吸收邊界反射，減少邊界反射對有效波場的干擾；疊前逆時偏移方法可有效識別碳酸鹽巖模型中的小尺寸孔洞和裂縫，歸為準(zhǔn)確，刻畫清晰，對碳酸鹽巖儲層成像具有顯著優(yōu)勢；GPU/CPU異構(gòu)并行加速技術(shù)可明顯提高偏移的計算效率，與傳統(tǒng)CPU算法相比，加速比可達60倍，對于實際推廣應(yīng)用具有重要意義。

碳酸鹽巖；疊前逆時偏移；PML；GPU/CPU異構(gòu)加速

碳酸鹽巖儲層多分布以裂縫及孔洞為油氣儲集空間。在我國已収現(xiàn)的油氣資源可采儲量及油氣產(chǎn)量中，碳酸鹽巖儲層只占可采儲量的18%，原油產(chǎn)量的3%，天然氣產(chǎn)量的27%，相較于全球的碳酸鹽巖油氣資源分布規(guī)律，我國碳酸鹽巖儲層的勘探和開収存在巨大的収展?jié)摿Α?/p>

本文采用基于聲波方程的疊前逆時偏移方法對碳酸鹽巖類構(gòu)造的儲層進行成像，作為目前最為精確的偏移算法，其采用雙程波動方程，能夠?qū)Σ▓鲋械亩喾N波（多次波、回轉(zhuǎn)波等）進行成像，并且對于復(fù)雜的地質(zhì)構(gòu)造具有更好適應(yīng)性。在偏移過程中，結(jié)合優(yōu)化的高階有限差分系數(shù)，在不增加計算量的前提下，提高數(shù)值模擬的精度，有效的壓制了數(shù)值頻散現(xiàn)象；引入完全匹配層（PML）邊界條件來吸收由計算機有限存儲所帶來的邊界反射，大大減少了邊界反射對有效波場的干擾；同時采用GPU/CPU異構(gòu)并行加速技術(shù)對算法進行優(yōu)化，在較大程度上提高了偏移的計算效率。通過對模型進行

測試，表明疊前逆時偏移可以對碳酸鹽巖類儲層中的孔洞及裂縫定位準(zhǔn)確，成像清晰，在成像上具有顯著優(yōu)勢，驗證了理論的準(zhǔn)確性和算法的可行性。

1 基本原理

1.1 疊前逆時偏移

疊前逆時偏移方法采用雙程波動方程，通過對波場中的多種波進行成像來達到對復(fù)雜構(gòu)造準(zhǔn)確成像的目的。其基本原理如下，根據(jù)實際數(shù)據(jù)確定觀測系統(tǒng)，然后依據(jù)得到的震源坐標(biāo)和檢波點坐標(biāo)分別模擬震源波場的正向傳播和檢波點處波場的反向傳播。在上述過程中，對每一時刻的震源波場和檢波點波場應(yīng)用成像條件進行成像，并對多炮結(jié)果進行疊加和去噪，最后獲得疊加剖面。

本文所采用的聲波方程如下所示：

式中U為相應(yīng)時刻的波場值，V為速度場，t為時間變量，x、z分別為水平和垂直方向上的空間變量。通過有限差分對其進行離散，離散后的形式［1］如下：

式中al為二階導(dǎo)數(shù)偶數(shù)階有限差分系數(shù)，具體推導(dǎo)公式如下：

差分系數(shù)對波場的數(shù)值模擬有著至關(guān)重要的影響，在這里，本文引入張金海［2］提出的優(yōu)化差分系數(shù)。與傳統(tǒng)差分系數(shù)相比，同階的優(yōu)化系數(shù)在計算量不變的前提下，可有效壓制由于數(shù)值離散引起的數(shù)值頻散，對波場模擬的精度更高。如圖1所示，上半部分為采用傳統(tǒng)差分系數(shù)的波場模擬結(jié)果，下半部分為采用優(yōu)化差分系數(shù)的結(jié)果，可以看出，相較于上半部分，下半部分的頻散現(xiàn)象得到了明顯的壓制，精度更高。

1.2 疊前逆時偏移

本文采用PML完全匹配層吸收邊界條件，能夠有效的壓制邊界反射，消除邊界干擾，基本實現(xiàn)原理是在傳統(tǒng)聲波方程中加入衰減算子，通過在有效波場的外部添加邊界衰減層來對邊界反射進行吸收衰減，衰減因子如圖2所示。當(dāng)波場在邊界內(nèi)部傳播時，衰減因子為零，方程變?yōu)樵镜男问剑瑹o吸收衰減，當(dāng)?shù)卣鸩▊髦吝吔鐓^(qū)域時，其能量將隨著衰減因子的變化而逐漸衰減，以達到吸收衰減邊界反射的目的。

引入衰減因子的聲波方程[5]如下：

式中d( x, z)為衰減因子，本文采用的衰減函數(shù)為余弦型衰減函數(shù)［3］：

式中B為衰減系數(shù)，PML為衰減層厚度，i為吸收點距離邊界處的距離。從式(4)可以看出，衰減函數(shù)有值時，該方程可以逐步的對波場能量進行衰減，當(dāng)衰減函數(shù)為零值時，方程退化為傳統(tǒng)的聲波方程，不影響波的正常傳播。波場衰減效果如圖3所示。

成像方面本文采用成像穩(wěn)定的互相關(guān)成像條件，將激収點的正傳波場與檢波點的反傳波場進行同時刻的互相關(guān)運算，公式如下：

式中S( x, z, t)為激収點 t時刻正傳波場，R( x, z, t)為檢波點t時刻反傳波場，I( x, z)為成像剖面。

1.3 GPU/CPU異構(gòu)加速運算

傳統(tǒng)的圖形處理器(GPU)專職用來進行2D、3D圖形渲染計算。隨著計算機技術(shù)的収展，GPU的強大的計算能力已遠(yuǎn)不止僅限于圖形計算，其在通用計算方面的優(yōu)勢已經(jīng)成功的應(yīng)用到了多項領(lǐng)域［7-10］，包括地震疊前時間偏移、基于SRME方法的多次波預(yù)測技術(shù)等。針對逆時偏移在帶來精確地成像結(jié)果的同時也引入了巨大的計算量的問題，運用GPU/CPU異構(gòu)加速運算可以得到良好的解決。

本文采用 Nvidia公司推出的 CUDA (Compute Unified Device Architecture)通用并行計算架構(gòu)。程序執(zhí)行時，將需要進行大規(guī)模計算的數(shù)據(jù)傳輸至顯存端，然后根據(jù)GPU的眾核結(jié)構(gòu)，將GPU運算核心劃分為眾多線程單元，每個線程負(fù)責(zé)一個數(shù)據(jù)單元的運算，最后所有線程同時并行計算。下面給出一個組對比代碼示例，代碼實現(xiàn)的是簡單的數(shù)組求和運算。傳統(tǒng)的CPU版本與GPU版本程序?qū)Ρ瓤梢?/p>

看出，采用GPU/CPU異構(gòu)并行運算可以利用多線程并行運算取代傳統(tǒng)CPU端冗雜的循環(huán)運算，可以大幅度的提高計算效率。

CPU版本代碼：

void sample_cpu(float **a_host,float **b_host,float **c_host,int width,int height)

{for(i=0;i＜height;i++)

for(j=0;j＜width;j++)

{c_host[i][j]=a_host[i][j]*b_host[i][j]}

}

GPU版本代碼：

__global__ void sample_gpu(float **a_device,float **b_device,float **c_device,int width,int height)

{int ty=blockIdx.y*blockDim.y+threadIdx.y;

int tx=blockIdx.x*blockDim.x+threadIdx.x;

c_device[ty*width+tx]=a_device[ty*width+tx]*b_d evice[ty*width+tx];

}

2 模型試算

本文根據(jù)碳酸鹽巖儲層的地質(zhì)構(gòu)造特點，設(shè)計了碳酸鹽巖地質(zhì)模型并進行疊前逆時偏移測試。

圖4所示的是碳酸鹽巖單洞模型，模型共設(shè)200個激収點，200道接收，炮間距、道間距均為10 m，網(wǎng)栺間距為10 m。其中孔洞直徑分別為10、20、30 m。對于成像過程中出現(xiàn)的低頻噪音，采用了拉普拉斯算子去噪方法進行壓制，本文就不在累述，最終偏移去噪結(jié)果如圖5所示。

圖6所示的是碳酸鹽巖裂縫模型，基本參數(shù)與碳酸鹽巖單洞模型相同。其中裂縫的直徑分別為10、20、30 m；裂縫長度分別為50、100、150 m。圖7為該模型的偏移去噪結(jié)果。

參與本文測試的計算機主要配置參數(shù)為：

CPU:Intel i3 2120;

內(nèi)存:2*4G DDR3 1600;

硬盤:希捷1TB 7200 r/s;

顯卡（GPU）:Nvidia Geforce GTX560，1 024 M顯存，顯存位寬256 bit,核心頻率850 MHz，顯存頻率4 500 MHz，流處理器個數(shù)336個。

3 結(jié) 論

疊前逆時偏移方法可以對地下微幅及大角度構(gòu)造進行精確成像。在碳酸鹽巖單洞模型的測試中，疊前逆時偏移方法對小尺寸孔洞成像明顯，歸位準(zhǔn)確；在碳酸鹽巖裂縫模型測試中可以看出，疊前逆時偏移方法在對小尺寸裂縫進行準(zhǔn)確識別的同時，對裂縫大角度分布情況也可以清晰的刻畫出來。在碳酸鹽巖構(gòu)造成像方面，疊前逆時偏移方法有著巨大的優(yōu)勢。

在模擬計算過程中，本文采用的優(yōu)化差分系數(shù)及PML吸收邊界條件保證了理論模擬的精確性，有效的壓制消除了頻散及邊界反射等噪音干擾。

GPU/CPU異構(gòu)加速技術(shù)為疊前逆時偏移計算提供了較高的計算性能，解決了傳統(tǒng)疊前逆時偏移算法中計算時間成本大的問題。在實測過程中，GPU/CPU異構(gòu)加速運算方案較傳統(tǒng)CPU運算方案可以取得60倍左右的加速比。

針對大規(guī)模實際數(shù)據(jù)處理，本文方法可引至高性能計算集群上進行實際應(yīng)用，具有廣闊的収展前景和實際應(yīng)用價值。

圖1 差分系數(shù)對比效果Fig.1 Difference coefficient comparison

圖2 衰減系數(shù)分布Fig.2 Attenuation coefficient distribution

圖3 波場衰減效果Fig.3 Wavefield attenuation effect

圖4 碳酸鹽巖單洞模型Fig.4 Carbonate single hole model

圖5 碳酸鹽巖單洞疊前逆時偏移去噪結(jié)果Fig.5 Carbonate single hole prestack reverse-time migration denoising results

圖6 碳酸鹽巖裂縫模型Fig.6 Carbonate fracture model denoising result

圖7 碳酸鹽巖裂縫疊前逆時偏移去噪結(jié)果Fig.7 Carbonate fracture prestack reverse-time migration

［1］劉洋, 李承楚, 牟永光. 仸意偶數(shù)階精度有限差分?jǐn)?shù)值模擬[J]. 石油地球物理勘探, 1998, 33(1): 1-10.

［2］Jin-Hai Zhang, Zhen-Xing Yao. Optimized finite-difference operator for broadband seismic wave modeling[J]. GEOPHYSICS, 2013, 78(1): 511-526.

［3］陳可洋. 完全匹配層吸收邊界天劍研究[J]. 石油物探, 2010, 49(5): 472-477.

［4］劉紅偉, 劉洪, 鄒振, 崔永福. 地震疊前逆時偏移中的去噪與存儲[J].地球物理學(xué)報,2010,53(9):2171－2180.

［5］郭念民, 吳國忱. 基于PML邊界的變網(wǎng)栺高階有限差分聲波方程逆時偏移[J]. 石油地球物理勘探, 2012, 47(2): 256-265.

［6］胡昊, 劉伊克, 常旭,等. 逆時偏移計算中的邊界處理分析及應(yīng)用[J]. 地球物理學(xué)報, 2013, 56(6): 2033-2042.

［7］李博, 劉紅偉, 劉國峰,等. 地震疊前逆時偏移算法的CPU/GPU實施對策[J]. 地球物理學(xué)報, 2010, 53(12): 2938-2943.

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

［9］張舒, 褚艷丼 ,趙開勇, 張鈺勃. GPU高性能運算之CUDA[M].北京：中國水利水電出版社,2009:0－276.

［10］張清, 遲旭光, 謝海波, 趙開勇, 吳慶, 陳維, 王獅虎, 褚曉文. 基于GPU實現(xiàn)疊前時間偏移走時計算的并行算法[J].計算機系統(tǒng)應(yīng)用, 2011, 20(8): 42-46.

Carbonate Reservoir Pre-stack Reverse-time Migration and GPU/CPU Heterogeneous Computing Research

WANG Kuo，CHEN Hong-ling，GUO Cong，ZHANG Wen-wu，YE Fei
（Northeast Petroleum University, Heilongjiang Daqing 163318，China）

The prestack reverse-time migration based on acoustic wave equation was adopted in imaging carbonate reservoir. In the process of migration, optimized differential coefficient and perfectly matched layer (PML) boundary condition were used to improve the accuracy of numerical simulation, and GPU/CPU heterogeneous parallel acceleration technology was also used to improve calculation efficiency. Model tests indicate that, compared with the traditional differential coefficient, the optimized coefficient with the same order has higher accuracy,and dispersion is suppressed obviously; PML boundary condition can absorb boundary reflection effectively and reduce interference effectively. Prestack reverse-time migration can identify small carbonate holes and cracks in the model, depict them clearly and accurately, it has significant advantage in imaging carbonate reservoir;GPU/CPU heterogeneous parallel technique can obviously improve the migration computation efficiency, compared with the traditional CPU algorithm, speed ratio can reach 60 times, is of great significance for actual application.

Carbonate; Prestack reverse-time migration; PML; GPU/CPU heterogeneous accerelating

TE 122

A

1671-0460（2014）10-2042-04

黑龍江省創(chuàng)新創(chuàng)業(yè)訓(xùn)練項目，項目號：201310220019。

2014-09-13

王闊（1992-），男，黑龍江大慶人，2011年考入東北石油大學(xué)地球物理專業(yè)，研究方向：從事地震資料處理方面的工作。E-mail：542093479@qq.com。