基于照明補償?shù)淖儭白硬ā蹦M成像方法

焦峻峰，趙愛國，廉西猛，崔慶輝，孫成禹*，吳涵

（1.中國石油大學（華東）地球科學與技術(shù)學院，山東青島 266580； 2.中國石化勝利油田物探研究院，山東東營 257022； 3.中山大學地球科學與工程學院，廣東珠海 519080）

0 引言

隨著勘探對象日益復雜且成像精度要求不斷提高，疊前深度偏移已成為目前地震勘探數(shù)據(jù)處理的重要技術(shù)。偏移可使繞射波收斂、反射波歸位，得到真實地反映地下構(gòu)造的結(jié)果。相比于時間域，疊前深度偏移對于復雜構(gòu)造具有更強適應(yīng)性，但深度域模型對于速度模型的準確度依賴性較高且運算量巨大。近年來，隨著計算機技術(shù)的飛速發(fā)展，計算能力不足已不再是深度偏移技術(shù)的瓶頸，該方法得以推廣應(yīng)用。為了提高計算效率，Micikevicius［1］在2009 年利用GPU 實現(xiàn)了有限差分的快速算法；李博等［2］、劉守偉等［3］及Foltinek等［4］相繼利用GPU 加速的有限差分算法實現(xiàn)了疊前逆時偏移成像，取得了一些成果。盡管如此，傳統(tǒng)的疊前深度偏移方法仍然存在正演應(yīng)用頻率高且運行周期長、疊前炮記錄偏移流程長等諸多問題。即便是模擬研究，仍需先正演后偏移的一系列繁雜流程。因此，對相關(guān)的應(yīng)用和研究極為不利。

成像算子在光學領(lǐng)域被定義為成像系統(tǒng)對點光源的脈沖響應(yīng)。理想的光學系統(tǒng)是物體上一個點發(fā)出的光，不經(jīng)歷任何散射、折射等作用，直接成像在一個點上。關(guān)于成像算子，最早可追溯到20 世紀90 年代，即Gelius 等［5-7］提出的廣義散射層析方法。隨后，Hamran等［8］實現(xiàn)了基于局部平面波的散射層析成像，Sethian［9］提出基于波前擴展的快速推進法射線追蹤方法。至此，基于成像算子的成像方法已初具雛形。

進入21世紀，張雪建等［10］提出在深度域合成地震記錄的方法； Lecomet等［11-14］綜合前人研究成果，提出了模擬疊前局部成像方法； 張金陵等［15］利用地震子波取得一維點擴散函數(shù)； 吳涵等［16］提出基于成像算子的模擬疊前深度成像方法； 段偉國等［17］提出基于密度擾動的兩步正演法計算點擴散函數(shù)（光學上表示一點產(chǎn)生的光源）。至此，該方法已臻于成熟，但尚存一些不足。如基于成像算子的模擬成像方法目前主要采用分區(qū)重構(gòu)和近似代替方法，在子波形態(tài)、能量歸位等方面難如人意； 同時，不可避免地出現(xiàn)突變界面，這些問題亟待解決。

針對深度域子波及地震資料，張雪建等［10］討論了地震子波時深轉(zhuǎn)換問題； 何惺華［18］對深度域子波、褶積和傅里葉變換等基本問題進行了充分討論和分析；林柏香等［19］將時間域計算方法引入深度域，得到深度域合成地震記錄； 樊中海等［20］分類介紹了地震反演方法，闡述了子波提取的重要性；邵廣周等［21］使用加窗自適應(yīng)子波提取方法從觀測記錄中提取與實際數(shù)據(jù)匹配較好的震源子波。

在此基礎(chǔ)上，本文發(fā)展了一種基于照明補償?shù)淖儭白硬ā蹦M疊前深度成像方法。該方法僅需速度模型、觀測系統(tǒng)及震源即可得到成像結(jié)果，簡化了先正演后偏移的傳統(tǒng)流程，提高了計算效率，且能適用于地下復雜構(gòu)造，對于觀測系統(tǒng)優(yōu)化及傳統(tǒng)成像結(jié)果分析具有重要指導意義。

1 模擬成像方法原理

常速介質(zhì)中，地震波頻率與波數(shù)的對應(yīng)關(guān)系為

式中：f為子波主頻；k為波數(shù)；v為介質(zhì)速度。

一維情況下，點擴散函數(shù)的波長與周期滿足

式中：T為時域子波的周期；λ為深度域點擴散函數(shù)的波長；v（x）為x處地層的速度。

給定時間域子波和速度，可得深度域點擴散函數(shù)（圖1a）； 將深度域點擴散函數(shù)變換到波數(shù)域，即可得波數(shù)域點擴散函數(shù)（圖1b）。

圖1 深度域點擴散函數(shù)（a）和波數(shù)域擴散函數(shù)（b）

分別用正演算子F和反演算子I表示地下介質(zhì)速度模型V的正演過程和成像過程，則正演結(jié)果和成像結(jié)果分別為

式中：F（r）為正演算子；I（r）為成像算子；S（r）為正演結(jié)果；R（r）為成像結(jié)果；V（r）為速度場；r為地下模型散射點位置。

如果將二者結(jié)合，最終成像結(jié)果可表示為

式中D（r）為直接成像算子。

已知標量聲波方程

式中： ?2為拉普拉斯算子；P(r，t)為地震波場。通過傅里葉變換，將式（5）變換到頻率域，同時引用目標函數(shù)，通過Born 近似可得到

式中：O（r）為地下速度擾動（可近似看作反射率）；V0（r）為背景速度場；Psc為散射波場；Pin為入射波場；rg為檢波點位置；rs為炮點位置。由于該背景場是不均勻且無衰減的，所以其格林函數(shù)可由射線追蹤得到［9］，同時結(jié)合Stamnes 給出的格林函數(shù)近似解［22］，引入偏移算子l(rg，rs，r')，得到成像結(jié)果（r'為成像點）

式中：D(r，r')表示直接成像算子；γ表示振幅；；。

利用Beylkin［23-24］等提出的嚴格偽逆理論可得

式中K＝-kd＋ku，其中kd為下傳波數(shù)向量（沿震源點到成像點的射線路徑的波數(shù)向量），ku為上傳波數(shù)向量（沿成像點到檢波器的射線路徑的波數(shù)向量）。若將直接成像算子趨近于δ函數(shù)，成像結(jié)果為

此時，只需將地震子波作用于直接成像算子，就可得到帶有子波響應(yīng)的模擬疊前深度剖面

式中：s(ω)為子波頻譜；K為照明信息參數(shù)。K的主要作用是提供照明點反射面對應(yīng)的傾角，即在地下一點處若有一個傾角與K垂直的反射面（圖2），則可進行模擬成像。

圖2 散射向量示意圖

直接成像算子可在波數(shù)域由點擴散函數(shù)與上傳波數(shù)向量和下傳波數(shù)向量構(gòu)建，通過計算下傳波數(shù)向量與上傳波數(shù)向量在成像點處的差值，可得到模擬成像所需的波傳向量K，因此波傳向量的構(gòu)建首先需計算成像點到各檢波點及震源點的射線路徑。以Marmousi模型（圖3a）為例，計算成像算子位于（2300 m，1420 m）的走時場（圖3b），得到該點射線路徑（圖3c），進而得到該點處的波傳向量（圖3d），最終與該點的點擴散函數(shù)在波數(shù)域構(gòu)建成像算子（圖3e）。

圖3 Marmousi 模型的成像算子求取過程

將地下速度模型的速度擾動信息與成像算子相互作用，即可得到模擬成像結(jié)果，因此只需將成像算子與反射系數(shù)模型做褶積即可。而空間域的褶積運算相當于波數(shù)域的乘積運算，所以需將反射系數(shù)模型（圖4a）變換到波數(shù)域（圖4b），與成像算子（圖3e）相乘，即可得到波數(shù)域成像結(jié)果（圖4d）； 再將波數(shù)域成像結(jié)果進行傅里葉反變換，即可得到模擬成像結(jié)果（圖4c）。

圖4 Marmousi 模型空間域反射系數(shù)模型（a）及其成像結(jié)果（c）、波數(shù)域反射系數(shù)模型（b）及其成像結(jié)果（d）

圖4 所示的Marmousi 模型試算，共461×294 個網(wǎng)格，采用正方形網(wǎng)格，網(wǎng)格大小為10 m。雷克子波主頻為30 Hz，采用地面放炮地面接收，共20 炮，每炮230道接收，道間距20 m?？梢娫摲椒ㄟm用性良好，模擬成像結(jié)果能較準確地反映地下構(gòu)造信息。但該方法所得成像結(jié)果在子波形態(tài)及能量分配等方面尚有待后續(xù)改善。

2 模擬成像法存在問題及解決方案

2.1 存在的問題

相較于傳統(tǒng)成像方法，模擬成像方法雖具有計算速度快、成像清晰的優(yōu)勢，但該方法主要利用成像算子與反射系數(shù)模型作褶積得到成像結(jié)果，無法達到子波隨地層速度而變化的效果，而成像算子的構(gòu)建若僅采用一點的子波信息，當算子位置不同時，由于所用子波不同，成像結(jié)果也會有所偏差。

同時，成像算子的構(gòu)建離不開波傳向量的計算。波傳向量提供的主要信息是潛在照明點反射面對應(yīng)的傾角，只有當該位置存在垂直于波傳向量的反射面，才能照明出該反射面。當算子位置不同或觀測系統(tǒng)改變時，波傳向量改變會導致照明情況不同，模擬成像結(jié)果也會存在差異。圖5展示了成像算子位置不同時的射線路徑及模擬成像結(jié)果。

圖5 不同位置處的射線路徑（上）及其模擬成像結(jié)果（下）

當成像算子在（100 m，2840 m）時（圖5d），僅有左下方一小部分能夠成像；當成像算子位于（2310 m，2840 m）處時（圖5e），模擬成像結(jié)果大致反映了地下構(gòu)造的基本形態(tài)，但在中部高陡傾角構(gòu)造處出現(xiàn)許多虛假構(gòu)造；當成像算子位于（4610 m，2840 m）處時（圖5f），深部的低速區(qū)成像結(jié)果不理想，出現(xiàn)缺失。這是因為成像算子位置不同時，射線路徑不同，通過射線路徑所構(gòu)造的波傳向量無法將地下構(gòu)造面照出，也就無法成像。

2.2 基于可變子波和照明補償?shù)哪M成像方法

獲取準確的深度域點擴散函數(shù)是模擬疊前成像方法的基礎(chǔ)。一個波動過程，既可表示為時間的函數(shù)，即時間域子波，表示某一質(zhì)點的各時刻振幅形成的曲線；也可表示為空間的函數(shù)——深度域點擴散函數(shù)，即一系列質(zhì)點的振動形成的曲線。一維情況下，點擴散函數(shù)的波長與周期關(guān)系由式（2）給出，給定時域子波和速度，通過式（2）即可得到深度域點擴散函數(shù)，常速介質(zhì)中頻率與波數(shù)的關(guān)系由式（1）給出。選擇時域內(nèi)主頻為30 Hz的雷克子波，得到對應(yīng)的時間域子波（圖6a），假設(shè)介質(zhì)速度均勻，分別給定介質(zhì)速度為3000 m/s 和1000 m/s，可得到對應(yīng)速度為3000 m/s（圖6b）和1000 m/s（圖6c）的深度域點擴散函數(shù)。

圖6 時間域子波（a）、介質(zhì)速度為3000 m/s（b）及介質(zhì)速度為1000 m/s（c）的深度域點擴散函數(shù)

由圖6b 和圖6c可見，地下介質(zhì)速度越大，深度域點擴散函數(shù)的波形更“胖”一些； 地下介質(zhì)速度越小，波形更“瘦”一些。深度域點擴散函數(shù)與地下介質(zhì)有關(guān)，是關(guān)于速度的函數(shù)，其形態(tài)、波數(shù)和相位都隨著速度的變化而改變； 由于“線性時不變”條件在深度域中并不成立，即在不同介質(zhì)分界面處波形并不對稱，因此使用單一的深度域點擴散函數(shù)是不合適的。

考慮使用多個點擴散函數(shù)，就須改變褶積方式。時間域褶積模型的實質(zhì)是時間域子波與對應(yīng)點的反射系數(shù)相乘相加的過程。而在深度域中，由于點擴散函數(shù)是速度的函數(shù)，褶積的過程可轉(zhuǎn)換為一種非穩(wěn)態(tài)的運算過程

式中：s為深度域合成記錄；wi為地下不同深度點的點擴散函數(shù)；Ri=(vi+1-vi)/(vi+1+vi)為反射系數(shù)。

圖7 為采用主頻為30 Hz 的雷克子波，在Marmousi速度模型中抽取第80 道（圖7a），按照時域子波以不同速度求取深度域點擴散函數(shù)，并與深度域的反射系數(shù)進行非穩(wěn)態(tài)褶積，得到模擬深度域的合成記錄（圖7b）。對比不同深度處的波形，可見其間的差距十分明顯。

考慮到成像算子的構(gòu)建需要波數(shù)域點擴散函數(shù)與波傳向量，而使用一點的點擴散函數(shù)會導致成像結(jié)果不準確，故采用單位脈沖替代波數(shù)域點擴散函數(shù)與波傳向量，構(gòu)建成像算子； 同時，在深度域中，反射系數(shù)模型與深度域點擴散函數(shù)作非穩(wěn)態(tài)褶積，得到深度域的模擬地震記錄，將深度域模擬記錄轉(zhuǎn)化到波數(shù)域并施加成像算子的作用，最后變換到深度域即可得到成像結(jié)果。具體計算流程如圖8所示。

圖8 模擬成像算法流程圖

同時，由于構(gòu)建成像算子時并未使用該點的子波而是用單位脈沖做替代，那么可將多個點的模擬成像結(jié)果疊加以解決常規(guī)模擬成像結(jié)果中能量不匹配或個別區(qū)域無法照明的問題。首先，需在深度域進行非穩(wěn)態(tài)褶積，得到模擬的地震記錄，然后將深度域模擬地震記錄變換到波數(shù)域，將其與成像算子相乘即可得到波數(shù)域成像結(jié)果，將波數(shù)域成像結(jié)果反變換，可得到變“子波”深度域模擬成像結(jié)果。

又因為傅里葉變換遵循線性性質(zhì)，所以模擬成像結(jié)果疊加運算可以轉(zhuǎn)換為成像算子疊加運算，也就是在波數(shù)域?qū)⒊上袼阕盈B加，然后將疊加后的成像算子作用在波數(shù)域模擬地震記錄上，即可得到波數(shù)域成像結(jié)果，進而將波數(shù)域成像結(jié)果反變換回深度域，得到最終的模擬成像結(jié)果。整個過程可表述為

式中：Di(r)為不同位置的成像算子；Ri(r)為不同位置處成像算子得到的模擬成像結(jié)果。

3 算例分析

按照本文提出的方法，采用主頻為30 Hz 的雷克子波，地面放30炮，共230道接收，以道間距20 m 進行計算，選用（100 m，2840 m）、（2310 m，2840 m）和（2310 m，2840 m）三點的成像算子疊加。得到照明補償后的模擬成像結(jié)果（圖9b）。與照明補償前模擬成像（圖9a）結(jié)果對比，兩種方法都可以較為準確地刻畫地下構(gòu)造形態(tài)，只有在一些細微的地方存在差異。

圖9 照明補償前（a）、后（b）的模擬成像結(jié)果

抽取單一成像算子模擬成像結(jié)果第80 道與深度域非穩(wěn)態(tài)褶積波形對比（圖10a），同時抽取照明補償后的變“子波”模擬成像結(jié)果第80 道波形與深度域非穩(wěn)態(tài)褶積波形對比（圖10b），可見改進后的模擬成像方法得到的成像振幅更接近理論振幅，成像結(jié)果中也包含了更多的細節(jié)，該方法不論在淺層還是深層成像振幅都與理論振幅吻合較好，相比單一成像算子模擬成像效果具有顯著提升。

圖10 使用單一算子（a）、疊加算子（b）的Marmousi 模型模擬成像結(jié)果波形與非穩(wěn)態(tài)褶積結(jié)果對比

在相同觀測參數(shù)下，采用不同方法分別得到Marmousi模型成像結(jié)果（圖11）。逆時偏移成像結(jié)果中振幅隨著深度增加顯著降低，導致深層一些構(gòu)造無法顯示，Kirchhoff偏移在復雜構(gòu)造或深部地區(qū)成像結(jié)果較差，對比發(fā)現(xiàn)改進的模擬成像方法能真實地反映地下構(gòu)造，為不同方法成像結(jié)果的評判提供參考。

圖11 不同成像方法結(jié)果對比

由于成像算子的構(gòu)建離不開射線路徑計算，所以，模擬成像結(jié)果的好壞可以間接反映觀測系統(tǒng)的選擇是否恰當。圖12 為地面230 道接收、道間距為20 m、炮數(shù)不同時的模擬成像結(jié)果。

圖12 Marmousi 模型不同炮點距時的模擬成像結(jié)果

對比圖12，地面啟爆30～40 炮，即炮間距為100～150 m 時模擬成像結(jié)果良好，模擬成像結(jié)果可作為傳統(tǒng)成像結(jié)果的參考，同時對于現(xiàn)場觀測系統(tǒng)的優(yōu)化具有指導意義。

4 結(jié)束語

本文改進了一種基于照明補償?shù)淖儭白硬ā蹦M疊前深度偏移成像方法。該方法可在已知速度模型的基礎(chǔ)下，根據(jù)觀測系統(tǒng)的布置，直接得到深度域模擬成像結(jié)果；相較于傳統(tǒng)方式大大節(jié)省計算時間，計算效率高； 同時，它在成像精度、振幅保持及子波形態(tài)方面有明顯優(yōu)勢，能精細地刻畫地下構(gòu)造形態(tài)，成像結(jié)果波形符合實際情況。因此，該方法對于指導觀測系統(tǒng)的優(yōu)化、成像結(jié)果評判及儲層預測等方面，均具有重要意義。