田坤, 王常波, 劉立彬,張建中

(1.中國(guó)石化勝利油田分公司 物探研究院, 山東 東營(yíng) 257022；2.中國(guó)海洋大學(xué) 海洋地球科學(xué)學(xué)院，山東 青島 266000)

0 引言

斜率層析成像是一種聯(lián)合使用反射波局部相關(guān)同相軸的走時(shí)和斜率數(shù)據(jù)，同時(shí)反演地層速度和反射點(diǎn)位置的有效方法[1]。該方法在共炮道集和共檢波點(diǎn)道集拾取反射波走時(shí)和斜率數(shù)據(jù)時(shí)，在層析反演中引入地震斜率信息，能夠克服反射波射線的多路徑問題，特別是在層析反演中不需要建立反射同相軸與地層界面的對(duì)應(yīng)關(guān)系，對(duì)低信噪比資料有更強(qiáng)的適應(yīng)性。斜率層析成像方法由Billette等[1]于1998年提出后的幾年內(nèi)，成功應(yīng)用于海上實(shí)際資料處理，并將二維斜率層析成像發(fā)展至三維情形[2-5]。斜率層析成像方法不斷得到完善和發(fā)展，提出了適于OBN資料反演的斜率層析成像[6]、各向異性斜率層析成像[7-8]、基于程函方程的伴隨狀態(tài)法斜率層析成像[9]、適用于起伏地表的斜率層析成像[10]、慢度平方三角網(wǎng)格斜率層析成像[11]等。在斜率層析成像的正則化方面，發(fā)展了基于反射角的結(jié)構(gòu)性平滑約束[12]、地層格架正則化約束[13]等方法。走時(shí)和斜率數(shù)據(jù)的拾取技術(shù)得到較大發(fā)展，一是在不同的數(shù)據(jù)域進(jìn)行拾取，如疊前時(shí)間域[14-15]、疊后時(shí)間域[16]以及深度偏移域[17-19]等數(shù)據(jù)域；二是對(duì)拾取的方法原理進(jìn)行改進(jìn)，如基于直線和雙曲線疊加的斜率自動(dòng)拾取[20]、基于曲波變換的地震反射波同相軸斜率的拾取[21]、應(yīng)用運(yùn)動(dòng)學(xué)反偏移提取斜率數(shù)據(jù)[22]等方法，提高了斜率數(shù)據(jù)的拾取精度。

對(duì)于實(shí)際地震資料，由于地震波的球面擴(kuò)散和介質(zhì)吸收衰減作用，使來自深部地層的地震波能量較弱[23-24]。在構(gòu)造復(fù)雜多樣、地層速度變化較大的地區(qū)，深部反射波變化復(fù)雜，地震資料信噪比往往較低。濟(jì)陽坳陷中—古生界潛山構(gòu)造復(fù)雜，疊前反射地震同相軸能量弱，十分雜亂，能夠拾取的深部反射波走時(shí)和斜率數(shù)據(jù)往往較少，導(dǎo)致深部地層射線覆蓋不足，層析反演速度的效果較差。為了解決該問題，本文在層析反演方程中，對(duì)數(shù)據(jù)關(guān)于地層速度的核函數(shù)進(jìn)行深度加權(quán)處理，增強(qiáng)深部反射波數(shù)據(jù)對(duì)深層速度的約束作用，從而提出一種基于深度加權(quán)的三維地震斜率層析成像方法，并通過理論模型合成的地震數(shù)據(jù)和濟(jì)陽坳陷古潛山實(shí)測(cè)地震資料的測(cè)試，驗(yàn)證了提出方法的有效性。

1 方法原理

地震斜率是反射波走時(shí)關(guān)于水平距離的導(dǎo)數(shù)，即慢度矢量的水平分量，可用之確定射線從炮點(diǎn)出發(fā)和到達(dá)檢波點(diǎn)的角度。地震立體層析成像[1]把一個(gè)炮—檢對(duì)的反射射線看成是從反射點(diǎn)分別到炮點(diǎn)和檢波點(diǎn)的兩段射線(如圖1所示)，通過不斷修改速度模型和反射點(diǎn)位置，使計(jì)算的兩段射線與地面的交點(diǎn)及對(duì)應(yīng)的斜率和走時(shí)與觀測(cè)數(shù)據(jù)達(dá)到一致，從而形成了一種有效的斜率層析成像方法[25]。為了提高對(duì)深部地層速度的層析反演效果，在地震立體層析成像方法基礎(chǔ)上，我們提出了一種基于深度加權(quán)的三維斜率層析成像方法。

圖1 模型參數(shù)與觀測(cè)數(shù)據(jù)幾何關(guān)系

1.1 模型離散化及走時(shí)和斜率計(jì)算

采用規(guī)則網(wǎng)格對(duì)速度模型進(jìn)行離散化，已知網(wǎng)格節(jié)點(diǎn)上的速度值，模型空間中任意位置的速度值由網(wǎng)格節(jié)點(diǎn)速度值的三次B樣條函數(shù)表示。當(dāng)給定射線的出發(fā)位置和出射角度時(shí)，采用二階Runge-Kutta法求解下列程函方程[2]，獲得從反射點(diǎn)到地表位置的射線路徑及其走時(shí)和斜率數(shù)據(jù)：

,

(1)

式中：x=(x,y,z)表示位置矢量，p=(px,py,pz)表示慢度矢量，t是沿射線波傳播時(shí)間，v為傳播速度，?代表梯度算子。

1.2 數(shù)據(jù)空間和模型空間

三維斜率層析的數(shù)據(jù)空間d可表示為：

,

(2)

其中，(xs,ys,zs)為炮點(diǎn)S的坐標(biāo)，psx和psy分別為炮點(diǎn)處沿x方向和y方向的斜率，(xr,yr,zr)為檢波點(diǎn)R的坐標(biāo)，prx和pry分別為檢波點(diǎn)處沿x方向和y方向的斜率，tsr為炮點(diǎn)S和檢波點(diǎn)R之間的反射走時(shí)，N表示射線對(duì)的個(gè)數(shù)，下標(biāo)n表示第n個(gè)射線對(duì)。炮點(diǎn)、檢波點(diǎn)的位置和高程從野外觀測(cè)系統(tǒng)文件中獲得，炮點(diǎn)、檢波點(diǎn)的斜率和雙程走時(shí)數(shù)據(jù)，用傾斜疊加等方法[14,20,21]在疊前地震道集上拾取。

模型空間m由離散速度參數(shù)mv和射線段參數(shù)mray組成，可表示為：

m=(mv,mray)T

,

(3)

(4)

(5)

其中，(xc,yc,zc) 為反射點(diǎn)C的坐標(biāo)，θs,φs,θr和φr分別表示從反射點(diǎn)到炮點(diǎn)S和檢波點(diǎn)R的射線出射角，ts和tr分別表示從反射點(diǎn)到炮點(diǎn)和檢波點(diǎn)的走時(shí)，N是炮-檢對(duì)或反射點(diǎn)的總數(shù)，vm是離散速度模型第m個(gè)節(jié)點(diǎn)的速度值，M表示整個(gè)離散速度模型的節(jié)點(diǎn)數(shù)。

1.3 反演方程及求解

斜率層析反演的目標(biāo)函數(shù)為：

(6)

其中，第一項(xiàng)是數(shù)據(jù)殘差項(xiàng)，第二項(xiàng)是模型約束項(xiàng)，λ是阻尼系數(shù)，調(diào)節(jié)這兩項(xiàng)作用的相對(duì)大小。dobs和dcal分別是觀測(cè)數(shù)據(jù)和計(jì)算數(shù)據(jù)，m是模型，m0是先驗(yàn)?zāi)Ｐ?，L表示拉普拉斯算子，對(duì)模型進(jìn)行光滑約束，上標(biāo)T表示矩陣轉(zhuǎn)置。Cd和Cm分別是數(shù)據(jù)和模型的協(xié)方差矩陣，分別對(duì)數(shù)據(jù)項(xiàng)和模型參數(shù)項(xiàng)進(jìn)行加權(quán)處理，并分別對(duì)走時(shí)、斜率等不同量級(jí)的觀測(cè)數(shù)據(jù)以及速度、反射點(diǎn)位置等不同量級(jí)的模型參數(shù)進(jìn)行歸一化處理。

dcal(m)是一個(gè)非線性函數(shù)，在初始模型m0處，對(duì)該函數(shù)進(jìn)行泰勒級(jí)數(shù)展開，取至線性項(xiàng)得

dcal(m)=dcal(m0)+G(m-m0)

,

(7)

其中，G=(Gv,Gray)T，Gv為數(shù)據(jù)關(guān)于離散模型節(jié)點(diǎn)速度的核函數(shù)矩陣，Gray為數(shù)據(jù)關(guān)于射線參數(shù)(反射點(diǎn)坐標(biāo)、射線出射角、射線走時(shí))的核函數(shù)矩陣[24]。

為了求取目標(biāo)函數(shù) (6)的極小解，把式 (7)代入式 (6)后，令目標(biāo)函數(shù)對(duì)模型空間參數(shù)的導(dǎo)數(shù)等于0，則得到下列線性反演方程:

，

(8)

其中，Δd表示數(shù)據(jù)殘差向量，其元素是觀測(cè)數(shù)據(jù)與當(dāng)前反演模型計(jì)算的數(shù)據(jù)之差，Δmv和Δmray分別是射線段參數(shù)和離散速度值的修正量。

反演方程(8)中包含各個(gè)反射波數(shù)據(jù)對(duì)其射線穿過的不同深度速度參數(shù)的核函數(shù)，我們對(duì)速度參數(shù)對(duì)應(yīng)的核函數(shù)進(jìn)行深度加權(quán)處理，建立下列含有深度加權(quán)的反演方程：

，

(9)

其中，wv是對(duì)速度參數(shù)對(duì)應(yīng)的核函數(shù)矩陣的深度加權(quán)矩陣，計(jì)算式見下面的式(11)，Δnray和Δnv分別是與射線段參數(shù)和速度參數(shù)對(duì)應(yīng)的待反演的未知向量。應(yīng)用LSQR算法[26]求解該反演方程，得到Δnray和Δnv。該解是引入深度加權(quán)后的解，再利用下式把該解換算成模型參數(shù)的修正量：

,

(10)

1.4 深度加權(quán)系數(shù)

斜率層析反演的模型空間包含反射點(diǎn)深度，每次迭代后獲得了各個(gè)炮—檢對(duì)的反射波數(shù)據(jù)對(duì)應(yīng)的反射點(diǎn)深度。用各個(gè)反射波數(shù)據(jù)對(duì)應(yīng)的反射點(diǎn)深度，建立式(9)中的深度加權(quán)系數(shù)矩陣wv。該矩陣是對(duì)角矩陣，其元素應(yīng)隨著離散速度節(jié)點(diǎn)深度的增大而增大。經(jīng)研究試驗(yàn)，構(gòu)造的深度加權(quán)系數(shù)表達(dá)式為：

，

(11)

式中，wik表示第i個(gè)反射波數(shù)據(jù)關(guān)于第k個(gè)節(jié)點(diǎn)速度的核函數(shù)的深度加權(quán)系數(shù)，zi是第i個(gè)反射點(diǎn)深度，hk是第i個(gè)反射波射線路徑穿過的第k個(gè)單元節(jié)點(diǎn)的深度，Ni是第i個(gè)反射波射線穿過的單元節(jié)點(diǎn)總數(shù)。

，

(12)

2 理論模型試驗(yàn)

根據(jù)濟(jì)陽坳陷埕島地區(qū)古潛山構(gòu)造特征和有關(guān)地震、鉆井資料，建立一個(gè)三維古潛山理論速度模型，用于測(cè)試本文方法的效果。如圖2所示，模型長(zhǎng)17.5 km，寬2 km，深7.8 km，速度變化范圍為2.2～8.0 km/s。淺部速度界面幾乎水平，且速度差異??；潛山頂位于4 km深度以下，潛山速度界面起伏較大，且與上覆地層存在較大的速度差異；潛山內(nèi)部地層無反射界面信息。將反射點(diǎn)均勻分布于各個(gè)速度界面，在各個(gè)反射點(diǎn)以θ為15°、30°、45°和55°，φ為10°、20°、30°和35°的角度向地表方向發(fā)射地震射線，追蹤射線至地表為止，獲得51 216組觀測(cè)數(shù)據(jù)。其中，來自潛山頂面以下的深部反射波數(shù)據(jù)較少。

圖2 古潛山理論速度模型

取y=1 km的速度模型切片(圖3a),初始速度v0=(2.2+0.25z)km/s(圖3b)。分別用無深度加權(quán)和有深度加權(quán)(α=1)的三維斜率層析成像方法對(duì)上述模擬數(shù)據(jù)進(jìn)行反演。模型離散網(wǎng)格在x、y和z這3個(gè)方向上的尺寸皆為250 m，反演速度模型在y=1km處的速度切片如圖3c和3d所示，反演速度沿深度的變化曲線與理論模型速度曲線的對(duì)比見圖4。可以看出，淺部地層反演模型與理論模型吻合較好，且有無深度加權(quán)，對(duì)潛山頂面上部地層的反演結(jié)果影響不大。但是，當(dāng)無深度加權(quán)時(shí)，對(duì)深層潛山的反演速度比理論模型速度整體偏小，且誤差較大；當(dāng)用了深度加權(quán)后，對(duì)深部潛山速度的反演效果得到了明顯提高，與理論速度變化趨勢(shì)吻合得更好。說明了層析反演中深度加權(quán)的作用及其有效性。

a—理論模型；b—初始模型；c—無深度加權(quán)的反演模型；d—有深度加權(quán)的反演模型

a—x=6 km處速度曲線；b—x=12 km處速度曲線

3 實(shí)際資料應(yīng)用

本文方法應(yīng)用于勝利油田孤島地區(qū)的三維地震資料。炮線和檢波點(diǎn)線分別呈EW向和SN向布設(shè)，炮間距50 m，檢波點(diǎn)間距25 m。對(duì)地震數(shù)據(jù)進(jìn)行噪聲壓制、振幅補(bǔ)償、初至切除等預(yù)處理。預(yù)處理后的單炮地震記錄，如圖5所示，其中1.5 s以上的淺部反射同相軸密集分布，而深部反射同相軸稀少，能拾取到的有效走時(shí)和斜率數(shù)據(jù)較少。基于傾斜疊加原理[14]，分別在共炮域和共檢波點(diǎn)域拾取了局部相干同相軸的走時(shí)和斜率，一共獲得150 047組觀測(cè)數(shù)據(jù)，用該組數(shù)據(jù)反演研究區(qū)的三維速度模型。

圖5 預(yù)處理后單炮地震記錄

根據(jù)研究區(qū)地表速度及有關(guān)速度等先驗(yàn)信息，建立了一個(gè)速度隨深度線性遞增的初始速度模型v0=(2.2+0.6z)km/s。模型在x方向長(zhǎng)度為15 km，在y方向?qū)挾葹? km，z方向厚度為5 km。模型離散網(wǎng)格在x、y和z方向的大小分別為250 m、250 m和200 m，分別使用無深度加權(quán)斜率層析成像方法和有深度加權(quán)斜率層析成像方法進(jìn)行了反演，深度加權(quán)系數(shù)中的α取值為0.5，迭代20次后反演的三維速度模型如圖6所示，沿y=9 km測(cè)線的反演速度模型切片如圖7所示?？梢钥闯?，兩個(gè)模型的淺層速度相差較小，而中深部速度相差較大，有深度加權(quán)斜率層析成像方法反演的潛山速度值明顯增大，且整個(gè)反演模型的變化細(xì)節(jié)更明顯。

a—無深度加權(quán)；b—有深度加權(quán)

a—無深度加權(quán)；b—有深度加權(quán)

分別使用無深度加權(quán)和有深度加權(quán)的斜率層析成像反演的速度模型進(jìn)行疊前深度偏移處理，沿y=9 km測(cè)線的疊前深度偏移剖面如圖8所示。比較兩個(gè)剖面可以看出，在淺部?jī)蓚€(gè)模型的偏移效果相差較小，但在中深部，特別在潛山頂面及其以下，使用深度加權(quán)斜率層析成像速度模型的偏移成像效果更好?？梢姡猩疃燃訖?quán)的斜率層析成像提高了反演速度模型的效果。

a—無深度加權(quán)斜率層析反演模型；b—有深度加權(quán)斜率層析反演模型

4 結(jié)論

當(dāng)?shù)卣鹂碧劫Y料信噪比較低時(shí)，能夠獲取的深層反射波走時(shí)和斜率數(shù)據(jù)往往較少，常規(guī)斜率層析成像對(duì)深層速度的反演效果較差。為此，本文利用離散模型節(jié)點(diǎn)深度和反演過程獲得的各個(gè)炮-檢對(duì)的反射點(diǎn)深度，構(gòu)造了一種針對(duì)斜率層析反演核函數(shù)矩陣的深度加權(quán)方程，通過對(duì)斜率層析反演方程進(jìn)行自適應(yīng)深度加權(quán)處理，使各個(gè)反射波數(shù)據(jù)對(duì)地下速度模型的約束作用隨著深度逐漸增強(qiáng)，從而提出一種基于深度加權(quán)的三維地震斜率層析成像方法。理論模型和實(shí)際資料測(cè)試結(jié)果表明，該方法能夠在保持淺層速度反演精度的同時(shí)，提高對(duì)深層速度的反演效果，改善整體反演速度模型的質(zhì)量，為低信噪比地震資料速度建模提供一種有效技術(shù)。