劉國峰 李海波 劉 語 孟小紅

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

0 引言

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

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

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

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

1 屬性偏移輔助波動(dòng)方程炮域偏移輸出炮檢距道集的理論與實(shí)現(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]指出，通過兩次偏移計(jì)算，可估算出加權(quán)于數(shù)據(jù)上的定量參數(shù)，即描述成如下過程

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

(2)

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

(3)

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

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

(4)

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

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

(5)

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

(6)

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

(7)

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

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

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

(8)

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

(9)

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

(10)

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

(11)

1.2 逆時(shí)偏移地表炮檢距輸出計(jì)算量估算

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

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

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

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

圖1 逆時(shí)偏移中兩次偏移輸出地表炮檢距道集的計(jì)算流程

2 正演模型試驗(yàn)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖11 用本文方法輸出的圖10標(biāo)識(shí)位置的地表炮檢距道集

4 結(jié)束語

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