屬性偏移輔助實現(xiàn)炮檢距域共成像點道集計算

劉國峰 李海波 劉 語 孟小紅

(①中國地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院，北京 100083； ②山東省煤田地質(zhì)局物測隊，山東濟南，山東 250104)

0 引言

勘探地震數(shù)據(jù)的成像方法包含Kirchhoff法和波動方程法(單程波、逆時偏移等)兩大類。Kirchhoff方法以單地震道為計算輸入，成像結(jié)果可根據(jù)輸入數(shù)據(jù)的炮檢距值進行排列以獲得炮檢距域共成像點道集，該道集中各道的炮檢距值來自于輸入道的炮檢點間距。為區(qū)別于其他局部炮檢距方法，也被稱為“地表炮檢距道集”。該道集是偏移迭代中速度更新的紐帶，可為地震數(shù)據(jù)AVO分析、屬性提取等深度信息發(fā)掘提供基礎(chǔ)數(shù)據(jù)[1]。

波動方程方法考慮了多值射線路徑等復(fù)雜傳播，具有更高的成像精度[2-6]。但炮域計算的波動方法偏移只能輸出零炮檢距剖面，無法應(yīng)用地表炮檢距道集實施類似Kirchhoff偏移的后續(xù)計算，帶來應(yīng)用上的局限性。如針對山地資料起伏地表進行波動方程偏移時，無法應(yīng)用Kirchhoff偏移，只能用波動方程偏移輸出的地表炮檢距道集進行速度更新[7-9]。部分波動方程偏移計算中，將炮點和檢波點位置的間距作為偏移輸出的標簽，在最終偏移結(jié)果輸出中重排以獲得一個炮檢距道集。該炮檢距并非地表炮檢距，不能進行速度更新等迭代。另外，波動方程偏移輸出的角度道集深層角度窄小[10-11]，且已有的研究[12]表明，不同方法計算的角度道集中，相同速度誤差導(dǎo)致的成像結(jié)果中同一成像點同相軸的曲率不同，表現(xiàn)出應(yīng)用于速度更新的不穩(wěn)定性。因此，如何在目前應(yīng)用的炮域波動方程偏移計算中高效輸出地表炮檢距道集，是具有實際應(yīng)用意義的研究方向。

炮域波動方程偏移輸出地表炮檢距道集最直接的方法是,將檢波點數(shù)據(jù)按炮檢距大小分成若干份后分別偏移，用炮檢距均值標識不同結(jié)果，得到地表炮檢距道集。該方法需對每炮數(shù)據(jù)多次執(zhí)行偏移計算，計算強度大，難以大規(guī)模實施[13]。Zhao等[14]改進了上述算法，只計算少數(shù)炮檢距段，再通過插值補齊所需炮檢距位置數(shù)據(jù)。另一類波動方程地表炮檢距輸出方法被稱為波動方程Kirchhoff偏移，因其所用Born方程與Kirchhoff積分公式相似，故而得名。差異之處在于：計算中炮點和檢波點傳播格林函數(shù)是用波動方程而非射線追蹤法[15-17]。

屬性偏移是通過兩次偏移計算，將地震數(shù)據(jù)中的幾何、物理屬性，如走時、反射角、入射角、炮檢距等，從地震數(shù)據(jù)映射到偏移結(jié)果中，從這兩次成像結(jié)果的比值提取屬性。該方法由Bleistein[18]提出，在立體層析等計算中獲得應(yīng)用[19]。Giboli等[20]和Lemaistre等[21]討論了應(yīng)用此方法輸出炮檢距屬性的可能。本文以屬性偏移為基礎(chǔ)，研究了應(yīng)用該方法進行地表炮檢距域成像點道集計算； 并以逆時偏移為例，分析了其計算強度(Computing intensity)； 分別以2D和3D模型對算法進行了數(shù)值驗證，并將其應(yīng)用于實際金屬礦地震數(shù)據(jù)的真地表成像計算，獲得了其地表炮檢距道集。

1 屬性偏移輔助波動方程炮域偏移輸出炮檢距道集的理論與實現(xiàn)

1.1 基于屬性偏移的炮域地表炮檢距道集理論

應(yīng)用Born方程[22-23]表示疊前炮域深度偏移

I(x,s)=?G(x,t|s,0)G(x,-t|r,0)×

(1)

Bleistein[18]指出，通過兩次偏移計算，可估算出加權(quán)于數(shù)據(jù)上的定量參數(shù)，即描述成如下過程

Iw(x,s)=?G(x,t|s,0)G(x,-t|r,0)×

(2)

式中w(x|s,r)是與成像點處炮點、檢波點位置等有關(guān)的加權(quán)參數(shù)，它直接作用于偏移前數(shù)據(jù)，可以是入射角、反射角、炮點坐標或檢波點坐標等。該加權(quán)參數(shù)可通過式(1)與式(2)兩次偏移結(jié)果的比值估算得到

(3)

若加權(quán)系數(shù)w(x|s,r)為地表炮檢距，則可計算出成像點位置的地表炮檢距值，進而輸出地表炮檢距道集。

將式(1)～式(3)簡化至波動方程偏移的互相關(guān)成像條件中[16]，具體表示為

(4)

式中：u(x,s,t)是t時刻的炮點波場；v(x,s,r,t)是t時刻s處激發(fā)的檢波點波場。

將炮檢距h=s-r與炮域檢波點數(shù)據(jù)相乘，再對乘積做偏移，結(jié)果為

(5)

通過求式(5)與式(4)比值，獲得單炮s在成像點x處的炮檢距為

(6)

根據(jù)獲取的成像區(qū)各點炮檢距值，將單炮成像結(jié)果排放到相應(yīng)的炮檢距中。依次計算所有炮點記錄并疊加，便可獲得最終的地表炮檢距道集

(7)

式中δ(·)為脈沖函數(shù)。

針對三維情形，依據(jù)式(4)～式(6)分別計算橫向和縱向炮檢距，繼而獲得總的炮檢距值。

據(jù)式(5)推知橫向炮檢距加權(quán)數(shù)據(jù)偏移結(jié)果為

(8)

同理，縱向炮檢距加權(quán)數(shù)據(jù)偏移結(jié)果為

(9)

據(jù)式(6)可得兩個方向的炮檢距分別為

(10)

根據(jù)求取的兩個方向炮檢距屬性，重復(fù)計算各炮，即可獲得最終的成像結(jié)果

(11)

1.2 逆時偏移地表炮檢距輸出計算量估算

波動方程偏移中，尤以逆時偏移計算強度最大。隨著計算技術(shù)的發(fā)展，特別是圖形處理器(GPU)的出現(xiàn)，通過眾核并行加速，大幅度提高了逆時偏移的有限差分計算效率[24]，推動了其工業(yè)計算應(yīng)用。在GPU上進行逆時偏移計算時，為了突破CPU和GPU間帶寬瓶頸的限制，采用隨機邊界處理條件[25]，兩次計算炮點波場傳播，以計算換存儲，節(jié)省了數(shù)據(jù)傳輸耗時。逆時偏移在GPU上的計算過程為： 炮點波場正傳到邊界處被隨機化，進行反傳； 檢波點波場反傳； 二者在相同時刻相關(guān)成像；逐炮偏移疊加，獲得整體的成像疊加剖面。

炮域波動方程偏移輸出的是零炮檢距剖面，即疊加剖面。根據(jù)上述屬性偏移理論，需通過調(diào)制后數(shù)據(jù)的偏移結(jié)果與原始數(shù)據(jù)的偏移結(jié)果的比值計算成像點處的炮檢距值，進而抽取獲得成像后的地表炮檢距道集。通過分析可知，兩種偏移結(jié)果除輸入的檢波點數(shù)據(jù)不同外，其他成像要素均一致，因此可將整個計算過程整合到一起。

總結(jié)成像過程為：炮點波場正傳到邊界處并被隨機化，然后反傳；檢波點波場和調(diào)制后的檢波點波場分別反傳，炮點波場和兩個檢波點波場分別在相同時刻相關(guān)成像，輸出兩種偏移單炮結(jié)果；應(yīng)用兩個波場比值抽取炮檢距道集，逐炮計算疊加，獲得最終成像的地表炮檢距道集結(jié)果。

圖1是逆時偏移輸出地表炮檢距道集計算流程。顯然，去掉調(diào)制波場成像部分即為普通逆時偏移。因炮點波場正、反兩次傳播未變，只是檢波點正傳部分增加了調(diào)制波場，故其總計算量僅增加約三分之一。

圖1 逆時偏移中兩次偏移輸出地表炮檢距道集的計算流程

2 正演模型試驗

選用Marmousi模型進行方法流程試算，并與Kirchhoff計算的地表炮檢距道集進行對比分析。

首先對所有240炮記錄進行炮檢距調(diào)制計算；然后分別對原數(shù)據(jù)和調(diào)制數(shù)據(jù)應(yīng)用處理流程(圖1)進行深度偏移，每個單炮記錄偏移結(jié)果根據(jù)式(6)所求炮檢距值進行炮檢距映射計算，疊加所有炮即可獲得最終地表炮檢距道集。

圖2是原始單炮記錄和地表炮檢距調(diào)制后的單炮記錄，從圖2b可見調(diào)制后單炮中地震道能量隨著炮檢距的增加而增大。

對原始數(shù)據(jù)和調(diào)制后的數(shù)據(jù)分別采用相同的成像參數(shù)進行偏移。針對原始單炮偏移結(jié)果(圖3a)，求取兩次偏移結(jié)果的比值，獲得成像點對應(yīng)的地表炮檢距值(圖3b)。

將成像結(jié)果按炮檢距值分選到不同的炮檢距，并將所有炮集記錄依據(jù)上述流程計算、疊加、輸出，獲得最終的地表炮檢距成像道集。應(yīng)用Kirchhoff深度偏移對Marmousi模型進行成像計算，將其地表炮檢距道集與本文方法進行對比。道集中最小炮檢距為200m，最大炮檢距為2600m，炮檢距步長為100m。圖4和圖5分別是CMP150和CMP350處兩種方法的道集對比。從中可見兩種道集的總體振幅特征相似，但因波動方程偏移成像質(zhì)量更好，淺層無動校拉伸問題，故儲層更能聚焦成像，道集整體更為平直和光滑。

圖2 Marmousi模型單炮記錄(a)及用地表炮檢距調(diào)制結(jié)果(b)

圖3 原始單炮偏移結(jié)果(a)及求取的成像點對應(yīng)炮檢距值(b)

當速度模型準確時，偏移獲得的共成像點道集呈平直狀態(tài)； 當速度模型存在誤差，道集會不同程度地向下彎曲和向上翹起，這也是速度更新的基礎(chǔ)。將Marmousi模型速度值分別提高和降低5%，進行輸出地表炮檢距道集計算。圖6顯示提速5%時，同相軸呈欠偏移而向下彎曲； 圖7是降速5%時，同相軸呈過偏移而向上彎曲。

采用窄方位角SEG/EAGE模型作為3D測試模型，共計4800炮。三維地表炮檢距計算時，需分別針對兩個方向進行單炮數(shù)據(jù)調(diào)制并偏移計算，獲取兩個方向的炮檢距分量后，采用式(12)抽取炮檢距道集，最終將多炮疊加而獲得最終道集剖面。

圖8是單炮記錄中部分排列原始記錄和兩個方向調(diào)制后的結(jié)果及其偏移結(jié)果的深度切片。圖9是從計算結(jié)果中抽取的單條疊加剖面和單個炮檢距道集。道集中，最小炮檢距為100m，最大炮檢距為2600m，炮檢距步長為100m。

圖4 CMP150處Kirchhoff偏移(a)與逆時偏移(b)輸出的地表炮檢距道集對比

圖5 CMP350處Kirchhoff偏移(a)與逆時偏移(b)輸出的地表炮檢距道集對比

圖6 提速5%欠偏移時成像道集剖面(a)及局部放大結(jié)果(b)

圖7 降速5%過偏移時成像道集剖面(a)及局部放大結(jié)果(b)

圖8 三維地表炮檢距道集計算中部分單炮記錄、調(diào)制結(jié)果及最終成像的深度切片a)部分原始單炮記錄； (b)橫向炮檢距調(diào)制后的部分單炮記錄； (c)縱向炮檢距調(diào)制后的部分單炮記錄； (d)原始單炮偏移結(jié)果深度切片； (e)圖b偏移結(jié)果深度切片； (f)圖c偏移結(jié)果深度切片

圖9 3D鹽丘模型(a)、對應(yīng)偏移剖面(b)及其中黑線處炮檢距道集(c)

3 實際數(shù)據(jù)計算

閩西南是中國重要的多金屬礦集區(qū)，推覆體是其中一種重要的控礦構(gòu)造。尋找深部推覆體賦存狀態(tài)是深部礦產(chǎn)勘探的重要內(nèi)容。地震勘探因其勘探深度和分辨率優(yōu)勢被寄予厚望。2017年，在該區(qū)采集了一條長約18000m二維地震測試剖面，炮間距為60m，道間距為20m，滿覆蓋次數(shù)約為60。該測線也是該區(qū)首條應(yīng)用于金屬礦勘探的地震測線。

區(qū)內(nèi)地表結(jié)構(gòu)復(fù)雜，高差大，基巖出露且地表不均勻性強，單炮記錄信噪比低，成像計算難度大。通過系列研究和對比，波動方程真地表偏移獲得最優(yōu)成像結(jié)果[26-27]，其中偏移速度模型來自地表層析反演速度與深層均方根速度轉(zhuǎn)換的層速度的疊合。深度偏移對速度模型反應(yīng)敏感，往往需幾輪偏移和速度迭代計算。目前，迭代計算及速度準確性判斷需要地表炮檢距道集，而本文研究使這種真地表波動方程偏移速度更新及判斷成為可能。

圖10a為地表與深層拼接的深度域速度模型；圖10b為應(yīng)用該模型進行起伏地表偏移的剖面，從該剖面中選擇三處(圖10b豎線)分別輸出炮檢距道，其最小炮檢距為100m，最大炮檢距為3600m，炮檢距步長為100m。

從圖11所示的地表炮檢距道集可看出，不同成像位置處的道集初至?xí)r間不同，代表了不同地表高程變化。道集上不同程度存在速度優(yōu)化的空間，同時道集上進一步的濾波等處理對提高最終成像質(zhì)量有潛在可能，這些問題有待后續(xù)研究。

圖10 波動方程真地表偏移速度模型(a)和偏移剖面(b)黑豎線標示炮檢距道集的三個位置

圖11 用本文方法輸出的圖10標識位置的地表炮檢距道集

4 結(jié)束語

本文介紹了一種基于屬性偏移的兩次偏移計算輸出地表炮檢距道集策略，它具有算法簡便、運算效率高的特點，適合于大規(guī)模數(shù)據(jù)應(yīng)用。但該方法計算的穩(wěn)定性要求所計算成像點僅有一個穩(wěn)相點，因此在復(fù)雜模型的炮檢距映射中需引入光滑或正則化計算等。波動方程地表炮檢距道集的高效計算，為單獨應(yīng)用波動方程偏移結(jié)果進行速度迭代和成像提供了基礎(chǔ)，為類似波動方程真地表偏移等特殊計算的速度迭代提供了保障。