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

劉國峰 李海波 劉 語 孟小紅

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

0 引言

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

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

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

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

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

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

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

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

(1)

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

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

(2)

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

(3)

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

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

(4)

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

將炮檢距h=s-r與炮域檢波點數據相乘，再對乘積做偏移，結果為

(5)

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

(6)

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

(7)

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

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

據式(5)推知橫向炮檢距加權數據偏移結果為

(8)

同理，縱向炮檢距加權數據偏移結果為

(9)

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

(10)

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

(11)

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

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

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

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

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

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

2 正演模型試驗

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 實際數據計算

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

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

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

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

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

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

4 結束語

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