基于GPU\\CPU協(xié)同并行計(jì)算的逆時(shí)偏移系統(tǒng)的實(shí)現(xiàn)策略

摘 要：在石油天然氣地震勘探數(shù)據(jù)處理方法中，逆時(shí)深度偏移以其近似條件要求低、成像精度高而一直是地震成像技術(shù)的發(fā)展方向。隨著PC—Cluster 集群的普及，逆時(shí)深度偏移的計(jì)算瓶頸開始得到解決，跨進(jìn)了真正工業(yè)應(yīng)用的門檻；而近年來通用用途GPU并行計(jì)算技術(shù)的發(fā)展，又為逆時(shí)深度偏移生產(chǎn)處理開辟了一個(gè)大幅改善性價(jià)比的前景。本文通過對(duì)逆時(shí)深度偏移成像條件，波場(chǎng)重構(gòu)的逆時(shí)遞推關(guān)系，以及GPU＼CPU協(xié)同并行計(jì)算方案的討論，介紹了基于GPU＼CPU協(xié)同并行計(jì)算的逆時(shí)偏移系統(tǒng)的實(shí)現(xiàn)策略。

關(guān)鍵詞：逆時(shí)深度偏析 ；GPU＼CPU協(xié)同；系統(tǒng)實(shí)現(xiàn)

1.前言

自從反射地震勘探進(jìn)入工業(yè)應(yīng)用以來，地震偏移一直是地震勘探資料成像處理方面最為重要的內(nèi)容。好的偏移處理技術(shù)對(duì)成像而言具有一錘定音的作用，因此地震偏移成像一直是地震處理方法研究中的一個(gè)熱點(diǎn)。新技術(shù)不斷推出，為日益復(fù)雜的地震勘探目標(biāo)提供了越來越精細(xì)、準(zhǔn)確的地震資料圖像描述。隨著計(jì)算機(jī)技術(shù)的高速發(fā)展，波動(dòng)方程偏移開始進(jìn)入疊前深度逆時(shí)偏移（RTM）的階段。

相對(duì)其他方法而言，逆時(shí)偏移有明顯的優(yōu)點(diǎn)。第一，它沒有傾角限制，可完成陡傾角反射層成像，而且在原理上可以利用轉(zhuǎn)換波、棱鏡波或多次反射波成像，并獲得更精確的振幅等動(dòng)力學(xué)信息，實(shí)現(xiàn)保幅成像。第二，其成像方法不受介質(zhì)速度變化的影響，可以更好地對(duì)復(fù)雜速度場(chǎng)進(jìn)行精細(xì)估計(jì)，能夠?qū)?fù)雜區(qū)域進(jìn)行較準(zhǔn)確的成像。

然而，逆時(shí)偏移技術(shù)也存在技術(shù)瓶頸，為了取得成像方面的普適性以及精細(xì)刻畫方面的優(yōu)勢(shì)，它要求的計(jì)算資源非常大。一方面，目前地球物理行業(yè)常規(guī)的PC‐Cluster 集群很難負(fù)荷承擔(dān)，另一方面，在價(jià)格成本方面付出的代價(jià)也非常高昂。為此，高效低成本的疊前深度逆時(shí)偏移整體方案仍然是目前技術(shù)研究上的一項(xiàng)重要內(nèi)容。而通用圖形處理芯片（GPU）計(jì)算技術(shù)的應(yīng)用為疊前深度逆時(shí)偏移處理技術(shù)提供了一種性價(jià)比較高的解決方案[1][2][3][4]。

2.基本原理

疊前逆時(shí)深度偏移利用接收到的地震記錄進(jìn)行波場(chǎng)重構(gòu)，并在重構(gòu)的波場(chǎng)中確定散射源的空間位置，利用散射源的空間分布構(gòu)筑波場(chǎng)的散射圖像，從而達(dá)到波場(chǎng)介質(zhì)成像的目的。換句話說，疊前逆時(shí)深度偏移剖面的本質(zhì)就是一幅波場(chǎng)散射源的空間分布圖。

根據(jù)惠更斯原理，人工激發(fā)的地震波傳播到介質(zhì)非均勻奇異點(diǎn)時(shí)，將發(fā)生散射現(xiàn)象。如果這種散射點(diǎn)成面狀（在二維空間中體現(xiàn)為線狀）連續(xù)分布時(shí)，在這個(gè)面上的各個(gè)散射點(diǎn)所生成的散射波將會(huì)疊加形成界面反射波。換句話說，反射是散射的一種特殊情況，而一切反射均可以經(jīng)由散射現(xiàn)象來描述。

因此，只要通過將各種散射波的能量歸位到各自散射源的空間位置上，就可以構(gòu)筑出介質(zhì)散射點(diǎn)的空間分布圖，從而，逆時(shí)深度偏移問題可以分解成如下兩個(gè)問題：

（1）波場(chǎng)散射能量的逆時(shí)反傳播問題；

（2）散射能量回傳到散射源時(shí)刻的確定（從而對(duì)波場(chǎng)進(jìn)行取樣成像）問題。

第一個(gè)問題可以通過波場(chǎng)正演模擬的逆過程來實(shí)現(xiàn)。對(duì)于給定的波場(chǎng)介質(zhì)空間，通過對(duì)介質(zhì)空間模型以及波場(chǎng)進(jìn)行離散化，可以用數(shù)字方法模擬波場(chǎng)的發(fā)生，傳播，一直到被記錄下來的全過程（詳見后文介紹）。利用這個(gè)模擬技術(shù)，也可以將在地表接收到的地震記錄作為震源信號(hào)對(duì)待，將其記錄過程逆過來，以結(jié)束時(shí)刻作為逆時(shí)傳播的開始時(shí)刻，模擬地震記錄重新被注入地下模型空間中并產(chǎn)生波場(chǎng)傳播的過程。這個(gè)過程中的直達(dá)波正是正過程散射波傳播的逆過程，在地震記錄上的散射波將在此過程中逐漸傳回原來的散射源所在位置。所以，這個(gè)過程是一個(gè)散射波場(chǎng)重構(gòu)的過程。

同時(shí)，除非極為簡(jiǎn)單的均勻介質(zhì)傳播情形并具備完美的邊界記錄，否則完整重構(gòu)波場(chǎng)的逆時(shí)傳播過程是不可能的。不過，就成像目標(biāo)而言，并不要求每個(gè)細(xì)微的波場(chǎng)現(xiàn)象都得到重構(gòu)，也不指望在重構(gòu)波場(chǎng)的過程中完全不產(chǎn)生新的人為干擾。只要這種被忽略掉的細(xì)微波場(chǎng)現(xiàn)象或產(chǎn)生的新的人為干擾不對(duì)最終的成像效果造成顯著影響，這種近似就是可以接受的。

第二個(gè)問題的解決在原理上相對(duì)簡(jiǎn)單，散射發(fā)生的時(shí)刻同時(shí)也是直達(dá)波到達(dá)的時(shí)刻，因此通過射線或者波動(dòng)方程模擬的方法確定直達(dá)波到達(dá)各個(gè)空間點(diǎn)的時(shí)刻，就可以用這個(gè)時(shí)刻對(duì)逆時(shí)重構(gòu)波場(chǎng)進(jìn)行采樣。如果該點(diǎn)曾經(jīng)發(fā)生過散射現(xiàn)象，則重構(gòu)波場(chǎng)在該點(diǎn)必有散射能量回歸，取樣的結(jié)果就能取到顯著的回歸到該點(diǎn)的散射能量。反之，取樣的結(jié)果將沒有能量，或者只取到很弱的隨機(jī)干擾的噪音的能量。最后，這樣取樣構(gòu)筑而成的剖面便可以反映出空間散射源的分布圖像，達(dá)到地震成像的目的。

采樣的方式有多種，最簡(jiǎn)單的方式是取成像時(shí)刻所對(duì)應(yīng)的波場(chǎng)值，但這存在一定的問題，在真實(shí)世界中，任何能量的發(fā)生過程都需要經(jīng)過一定的時(shí)間才能達(dá)到峰值。散射發(fā)生時(shí)刻的震源子波信號(hào)的能量還沒來得及積累起足夠的能量，這時(shí)候的子波幅值是比較小的，對(duì)這個(gè)幅值的采樣不能真實(shí)反映散射的發(fā)生。當(dāng)然，如果將采樣時(shí)間做適當(dāng)后延，使得采樣剛好取得子波能量最強(qiáng)的相位，這個(gè)問題好像是可以解決的，但是，地震子波在傳播過程中存在大地濾波作用，子波波形在傳播過程中會(huì)發(fā)生變化，因此無法在整個(gè)模型空間的采樣過程中采用一個(gè)統(tǒng)一的時(shí)延來調(diào)整采樣時(shí)間，使得采樣剛好可以采得子波的最大相位振幅。所以，這種用一個(gè)時(shí)刻的波場(chǎng)值來構(gòu)筑介質(zhì)圖像的方法在計(jì)算穩(wěn)定性方面是相當(dāng)脆弱的。

一個(gè)比較穩(wěn)定可靠的采樣方法是用發(fā)生在該空間點(diǎn)上的兩個(gè)波場(chǎng)的時(shí)間過程做互相關(guān)來實(shí)現(xiàn)，這在下面的成像條件一節(jié)詳細(xì)討論，并具體給出數(shù)學(xué)表達(dá)公式。而且我們將發(fā)現(xiàn)，上述用散射時(shí)刻的重建波場(chǎng)值來構(gòu)建散射圖像的方法，實(shí)質(zhì)上不過是互相關(guān)成像條件在相關(guān)時(shí)間序列長(zhǎng)度為零的特例而已。

3.實(shí)現(xiàn)策略

研究目標(biāo)是根據(jù)地震勘探觀測(cè)記錄得到的共炮點(diǎn)道集疊前數(shù)據(jù)，觀測(cè)系統(tǒng)以及波場(chǎng)傳播速度模型，通過逆時(shí)偏移計(jì)算，構(gòu)筑出地下介質(zhì)的空間分布圖像，這個(gè)目標(biāo)可以通過如下過程達(dá)到。

對(duì)于每一個(gè)地震炮集疊前數(shù)據(jù)，觀測(cè)系統(tǒng)及其波場(chǎng)傳播速度的深度模型，先在地震數(shù)據(jù)上提取根震源信號(hào)，并通過正演模擬計(jì)算合成參考波場(chǎng)的邊界時(shí)間記錄。

（1）在相同的速度模型上，用相同的逆時(shí)計(jì)算方式，平行同步地分別用觀測(cè)地震記錄重構(gòu)觀測(cè)波場(chǎng)和合成參考波場(chǎng)的邊界時(shí)間記錄重構(gòu)參考波場(chǎng)；

（2）在逆時(shí)遞推過程中，同步計(jì)算重構(gòu)的觀測(cè)波場(chǎng)與重構(gòu)的參考波場(chǎng)的互相關(guān)函數(shù)，作為散射信號(hào)源點(diǎn)的圖像描述；

（3）用上述方法對(duì)每炮資料逐一處理，并在剖面模型上進(jìn)行逐炮相關(guān)數(shù)據(jù)圖像的疊加，最終得到整條測(cè)線的疊前逆時(shí)深度偏移剖面。

上述方法流程具體如圖1所示。

即便采用二步法確定逆時(shí)參考波場(chǎng)可以大幅度節(jié)省計(jì)算資源，但疊前逆時(shí)深度偏移本身的計(jì)算量仍然非常巨大，因此可以采用GPU＼CPU協(xié)同的并行架構(gòu)來解決疊前逆時(shí)深度偏移的繁重計(jì)算任務(wù)。協(xié)同計(jì)算的分工可以設(shè)計(jì)為：主機(jī)CPU負(fù)責(zé)輸入輸出，頂層邏輯以及各種零碎的計(jì)算工作；計(jì)算量密集的波場(chǎng)遞推模擬計(jì)算過程由GPU負(fù)責(zé)，用CUDA編程實(shí)現(xiàn)。目前用于科學(xué)與工程方面進(jìn)行通用編程計(jì)算的NVIDIA的Tesla序列GPU產(chǎn)品基本上可以提供足夠的設(shè)備儲(chǔ)存空間進(jìn)行逐炮的二維逆時(shí)深度偏移計(jì)算。因此，在GPU上的逆時(shí)深度偏移并行計(jì)算主要在線程層面上執(zhí)行。

4.結(jié)論

本文介紹了基于GPU＼CPU協(xié)同并行計(jì)算的逆時(shí)偏移系統(tǒng)的實(shí)現(xiàn)策略，具有高效低成本的優(yōu)勢(shì)，主要在逆時(shí)深度偏移計(jì)算環(huán)節(jié)上，采用GPU＼CPU協(xié)同的并行架構(gòu)來解決疊前逆時(shí)深度偏移的繁重?cái)?shù)值計(jì)算任務(wù)，其中安排主機(jī)CPU負(fù)責(zé)輸入輸出，頂層邏輯以及各種零碎的計(jì)算工作，計(jì)算量密集的波場(chǎng)遞推模擬計(jì)算過程由GPU負(fù)責(zé)，用CUDA編程實(shí)現(xiàn)，GPU上的逆時(shí)深度偏移并行計(jì)算主要在線程層面上執(zhí)行。

參考文獻(xiàn)：

[1]付小波， 馬中高， 余嘉順， 原健龍， 韓超。基于多圖形處理單元加速的各向異性彈性波正演模擬[J].科學(xué)技術(shù)與工程，2018，18（11）：16-22.

[2]付小波，余嘉順，原健龍。各向異性彈性波正演的多GPU加速計(jì)算策略[C]. 2017中國(guó)地球科學(xué)聯(lián)合學(xué)術(shù)年會(huì)論文集（二十五）——專題50：地震波傳播與成像.

[3]張猛。基于GPU并行加速的黏聲最小二乘逆時(shí)偏移及應(yīng)用[J]. 巖性油氣藏，2022，34（01）：148-153.

[4]劉聰?；诙郍PU的逆時(shí)偏移并行化方法研究[D].湖南大學(xué)，2021.