李海山　楊午陽　雍學(xué)善

(①中國石油勘探開發(fā)研究院西北分院，甘肅蘭州 730020； ②CNPC油藏描述重點(diǎn)實(shí)驗(yàn)室，甘肅蘭州 730020)

1　引言

基于模型的波阻抗反演方法由于具有分辨率高、抗干擾性強(qiáng)等優(yōu)點(diǎn)，在相關(guān)領(lǐng)域得到廣泛應(yīng)用[1-6]。盡管如此，這種方法的反演結(jié)果嚴(yán)重依賴于初始模型[7]，因此初始模型的建立至關(guān)重要。目前工業(yè)界普遍采用的初始模型建立方法是以地震解釋層位作為橫向控制、以斷層和測井資料等作為縱向控制，根據(jù)設(shè)定的沉積模式將已知屬性進(jìn)行外推內(nèi)插[8]，這類方法的顯著缺點(diǎn)是需要精細(xì)的層位和斷層解釋成果，但在地質(zhì)構(gòu)造及沉積關(guān)系復(fù)雜的區(qū)域，往往得不到精細(xì)的層位和斷層解釋結(jié)果，且人工解釋不但耗時(shí)還很容易產(chǎn)生解釋誤差[9]。

為了克服層位約束建模法的缺陷，有人提出了數(shù)據(jù)驅(qū)動(dòng)的初始模型建立方法[10-18]，這類方法既不需要耗時(shí)的精細(xì)層位和斷層解釋，又充分利用了地震數(shù)據(jù)本身所隱含的構(gòu)造和沉積信息。Schultz等[10，11]、Ronen等[12]和Hampson等[13]通過利用地質(zhì)統(tǒng)計(jì)學(xué)、多元線性回歸或人工神經(jīng)網(wǎng)絡(luò)等建立的地震屬性與測井信息間的關(guān)系獲得波阻抗、孔隙度、泥質(zhì)含量等各種參數(shù)模型； Hale[14， 15]提出成像引導(dǎo)的建模方法，通過求解程函方程和混合鄰近插值方程得到符合地質(zhì)規(guī)律的參數(shù)模型； 劉賢紅等[16]提出基于地震屬性約束的初始波阻抗模型建立方法，在縱向上利用地震波形相關(guān)分析尋找最優(yōu)插值點(diǎn)，在橫向上將地震資料的振幅、波形相關(guān)系數(shù)等地震屬性轉(zhuǎn)換為權(quán)系數(shù)進(jìn)行插值建模； Fomel[17]和Karimi等[ 18]提出基于地震數(shù)據(jù)的預(yù)測充填建模方法，該方法首先利用平面波解構(gòu)法(plane-wave destruction，PWD)[19， 20]計(jì)算得到分別沿著縱測線和聯(lián)絡(luò)測線方向的同相軸局部斜率，然后根據(jù)同相軸局部斜率將井點(diǎn)的信息延拓到整個(gè)三維空間。

平面波解構(gòu)最初是由Claerbout[21]在1992年提出的，由于基于平面波解構(gòu)的方法具有簡單快速、計(jì)算精度高、穩(wěn)定可靠的優(yōu)點(diǎn)，目前已經(jīng)廣泛地應(yīng)用于地震資料處理、成像、解釋、反演等各個(gè)方面[22-24]。本文受Fomel[17]和Karimi等[ 18]研究工作的啟發(fā)，在對(duì)平面波解構(gòu)進(jìn)行深入研究的基礎(chǔ)上，提出了基于平面波解構(gòu)的初始波阻抗建模方法。首先利用平面波解構(gòu)法計(jì)算得到地層局部斜率，接著利用平面波解構(gòu)方程所隱含的道與道之間的預(yù)測關(guān)系，從井點(diǎn)位置開始進(jìn)行波阻抗信息的延拓?？紤]到直接進(jìn)行計(jì)算的不適定性，在具體計(jì)算時(shí)，引入Tikhonov正則化方法提高延拓算法的穩(wěn)定性及抗噪性。模型數(shù)據(jù)測試及實(shí)際資料應(yīng)用證明了本文方法的有效性，利用本文方法建立的初始波阻抗模型與地震記錄所反映的地質(zhì)規(guī)律具有較好的一致性，克服了層位約束建模法容易產(chǎn)生的局部地質(zhì)現(xiàn)象與解釋層位不匹配的缺陷。

2　基于PWD的局部斜率計(jì)算原理

2.1　平面波解構(gòu)算子

Claerbout于1992年提出了平面波解構(gòu)的概念，其主要思想來源于描述地震數(shù)據(jù)的局部平面波模型[23]。局部平面波微分方程可表示為[19]

(1)

式中：u=u(t,x)為波場；t和x分別為時(shí)間和道間距；σ=σ(t,x)為同相軸局部斜率。在斜率為常數(shù)的情況下，式(1)的通解可表示為[19]

u(t,x)=f(t-σx)

(2)

其中f(t)代表波場。若斜率σ不依賴于時(shí)間t，可將式(2)變換到頻率域，其微分方程形式為[19]

(3)

式(3)的通解為[19]

(4)

(5)

(6)

C(zt,zx)=A(zt,zx)B(1/zt)

=B(1/zt)-zxB(zt)

(7)

式中濾波器B(zt)的系數(shù)可以通過在零頻率附近進(jìn)行泰勒級(jí)數(shù)展開確定，其中通常采用的三點(diǎn)中心濾波器B3(zt)為[19]

(8)

式中的系數(shù)分別為

(9)

2.2　局部斜率計(jì)算方法

平面波解構(gòu)算子是一種時(shí)空域預(yù)測誤差濾波器，它通過利用相鄰地震道預(yù)測當(dāng)前地震道的方法估算同相軸局部斜率，其預(yù)測準(zhǔn)則是原始地震道與預(yù)測出的地震道之間的殘差能量最小。如果用C(σ)表示濾波器C(zt,zx)在時(shí)間域與地震數(shù)據(jù)u(t,x)進(jìn)行二維褶積運(yùn)算的算子，則局部平面波微分方程式(1)的隱式有限差分形式為

C(σ)u≈0

(10)

在最小平方意義下，局部斜率σ的估算可轉(zhuǎn)化為求解最小二乘目標(biāo)函數(shù)[19]

(11)

由式(9)可知最優(yōu)化問題(式(11))是局部斜率σ的非線性優(yōu)化問題，用高斯—牛頓法(式(11))可轉(zhuǎn)化為[19]

C′(σ0)Δσu+C(σ0)u≈0

(12)

式中： Δσ為斜率增量；σ0為初始斜率；C′(σ)為C(σ)對(duì)σ的偏導(dǎo)數(shù)。求解式(12)后，初始斜率σ0通過加上斜率增量Δσ得到更新； 然后循環(huán)迭代求解式(12)。該方法需要數(shù)次非線性迭代獲得期望的局部斜率估計(jì)，收斂速率依賴于給定的起始斜率值。為了提高計(jì)算方法對(duì)噪聲的抑制能力以及計(jì)算結(jié)果的連續(xù)性，在局部斜率更新過程中，可引入預(yù)條件方法[19， 22]

εDΔσ≈0

(13)

式中：ε為常系數(shù)擾動(dòng)因子；D為預(yù)條件算子。

3　基于PWD的初始波阻抗建模原理

3.1　初始波阻抗建模原理

計(jì)算得到同相軸局部斜率之后，就可以利用平面波解構(gòu)濾波器通過已知道的波阻抗值預(yù)測相鄰道的波阻抗值。如果采用式(9)給出的三點(diǎn)中心濾波器系數(shù)，則式(10)可寫成矩陣方程形式

C1ux-C2ux-1≈0

(14)

式中

(15)

(16)

ux=[u(1,x)u(2,x)u(3,x)…u(T,x)]T

(17)

ux-1=[u(1,x-1)u(2,x-1)u(3,x-1)…u(T,x-1)]T

(18)

式(15)～式(18)中T為最大的時(shí)間采樣點(diǎn)數(shù)(1≤t≤T)，系數(shù)

(19)

根據(jù)式(14)，如果由x道預(yù)測x-1道，則

(20)

如果由x道預(yù)測x+1道，則

(21)

根據(jù)式(20)和式(21)，如果利用鉆井所在位置x道的波阻抗值，就可計(jì)算得到相鄰的x-1道和x+1道的波阻抗值，依此類推，就可得到整個(gè)剖面的波阻抗值。如果有多口井資料，可利用如反距離加權(quán)等方式得到多口井約束的建模結(jié)果。

3.2　Tikhonov正則化

為方便表述，將式(20)和式(21)用如下的形式統(tǒng)一表示

Cu≈w

(22)

根據(jù)矩陣的奇異值分解理論，式(22)中的稀疏矩陣C可分解為列正交矩陣S、元素大于或等于零的對(duì)角矩陣Σ和正交矩陣V的轉(zhuǎn)置的乘積，即稀疏矩陣C可表示為

(23)

其中左奇異向量st和右奇異向量vt分別為矩陣S和V的正交列向量，奇異值滿足λ1≥λ2≥…≥λT。式(22)的最小二乘解可以表示為

(24)

由于稀疏矩陣C可能為病態(tài)矩陣，存在較小的接近于零的奇異值，用式(24)直接求解式(22)可能會(huì)使所求解遠(yuǎn)遠(yuǎn)偏離原始問題的真解。對(duì)于這種不適定問題的求解，需要采用有效的正則化方法。其中Tikhonov正則化法是最著名的一種[25，26]，其基本思想是將解的范數(shù)最小作為約束條件加入到原問題中，將原問題轉(zhuǎn)化為在數(shù)據(jù)擬合程度和解范數(shù)最小之間達(dá)到某種平衡的問題，該問題即如下最小化問題

(25)

式中μ為正則化參數(shù)。根據(jù)矩陣C的奇異值分解，正則化參數(shù)μ對(duì)應(yīng)的正則化解為

(26)

4　模型數(shù)據(jù)測試

為說明基于平面波解構(gòu)的初始波阻抗建模方法的性能，下面用模型數(shù)據(jù)進(jìn)行測試。圖1a為利用部分Marmousi 2模型得到的二維縱波合成地震記錄； 圖1b為利用平面波解構(gòu)法得到的局部斜率場，負(fù)值表示局部同相軸向左傾斜，正值表示局部同相軸向右傾斜(圖中轉(zhuǎn)換成了傾角，傾角θ與斜率σ的關(guān)系為θ=arc tanσ)。由于圖1a中的地震反射同相軸相對(duì)比較平緩，圖1b中大部分斜率值都接近于零，在局部位置出現(xiàn)了相對(duì)較大的斜率值，整體上與合成地震記錄對(duì)應(yīng)關(guān)系良好。

圖1c為采用Jason軟件建立的初始阻抗模型，建模時(shí)以由圖1a解釋得到的從淺到深的23個(gè)層位作為橫向控制，以圖中黑色曲線所示的縱波阻抗曲線作為縱向控制； 圖1d為采用本文方法建立的初始波阻抗模型，建模時(shí)采用圖1b所示的局部斜率場，以圖中黑色曲線所示的縱波阻抗曲線作為已知的道向兩邊的道進(jìn)行遞推計(jì)算。 由圖1c和圖1d可見，盡管本文方法建立的初始波阻抗模型在遠(yuǎn)離已知點(diǎn)處準(zhǔn)確性略有降低，但整體上與合成地震記錄所反映的地質(zhì)規(guī)律具有較好的一致性，且不需要像Jason軟件那樣費(fèi)時(shí)的精細(xì)層位解釋。

圖1　模型數(shù)據(jù)測試結(jié)果比較

圖1e為利用Jason軟件建立的初始阻抗模型反演得到的結(jié)果；圖1f為利用本文方法建立的初始波阻抗模型反演得到的結(jié)果。從反演結(jié)果可見，圖中黑色橢圓內(nèi)的低阻氣藏及黑色矩形內(nèi)的低阻油藏在兩種初始模型的反演結(jié)果上都得到了準(zhǔn)確的刻畫，且本文方法反演結(jié)果在局部細(xì)節(jié)上的刻畫能力優(yōu)于Jason軟件反演結(jié)果，這充分說明了本文方法的有效性和優(yōu)越性。

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

為進(jìn)一步說明本文初始波阻抗建模方法的性能，將本文方法應(yīng)用于實(shí)際資料的建模和反演中。圖2a為M油田的疊后地震記錄，在CDP314處有一口生產(chǎn)井，圖中的黑色曲線為該井的縱波阻抗曲線；圖2b為利用平面波解構(gòu)法得到的局部斜率場；圖2c為以圖2a中紅色線所示的兩個(gè)解釋層位作為橫向控制，以井點(diǎn)的縱波阻抗曲線作為縱向控制，采用與層位頂、底平行的沉積模式，由Jason軟件建立的初始阻抗模型；圖2d為采用本文方法建立的初始波阻抗模型；圖2e為利用Jason軟件建立的初始阻抗模型反演得到的結(jié)果；圖2f為利用本文方法建立的初始波阻抗模型反演得到的結(jié)果。

由圖2c可見，由于解釋層位較少，采用簡單的沉積模式很難反映目的層內(nèi)的真實(shí)地質(zhì)規(guī)律，因此用Jason軟件建立的初始波阻抗模型很容易出現(xiàn)穿時(shí)現(xiàn)象，初始模型不能客觀地反映地下地質(zhì)規(guī)律。而由圖2d可見，本文初始波阻抗建模方法由于是沿著局部同相軸傾斜方向延拓井點(diǎn)的波阻抗值，能夠得到與地震記錄所反映的地質(zhì)規(guī)律高度一致的結(jié)果。由圖2e和圖2f可見，盡管利用兩種不同的初始模型都準(zhǔn)確地反演出了圖中黑色橢圓內(nèi)的低阻氣藏，但在遠(yuǎn)離井點(diǎn)(紅色箭頭所指)位置，利用Jason軟件建立的初始阻抗模型反演得到的結(jié)果模型化現(xiàn)象嚴(yán)重，可靠性遠(yuǎn)地低于利用本文方法建立的初始波阻抗模型反演得到的結(jié)果。

圖2　實(shí)際資料應(yīng)用結(jié)果比較

6　結(jié)論

(1)本文首先利用平面波解構(gòu)法計(jì)算得到地層局部斜率，接著利用平面波解構(gòu)方程所隱含的道與道之間的預(yù)測關(guān)系進(jìn)行波阻抗信息的空間延拓，并引入Tikhonov正則化方法來提高延拓算法的穩(wěn)定性及抗噪性，實(shí)現(xiàn)了基于平面波解構(gòu)的初始波阻抗建模方法。

(2)本文方法充分利用了地震數(shù)據(jù)中所包含的構(gòu)造和沉積信息，不再像傳統(tǒng)建模方法那樣需要精細(xì)的層位及斷層解釋結(jié)果，直接利用地震數(shù)據(jù)和測井信息進(jìn)行初始模型的建立。模型數(shù)據(jù)測試及實(shí)際資料應(yīng)用證明了本文方法的有效性。

(3)本文方法建立的初始波阻抗模型與地震相變化具有較好的一致性，克服了由于地震層位解釋不準(zhǔn)確以及地層與解釋層位關(guān)系復(fù)雜造成的模型偏差現(xiàn)象，是對(duì)層位約束建模方法的一種補(bǔ)充和完善，通過優(yōu)化完善并推廣到三維初始模型的建立，具有廣闊的應(yīng)用前景。

[1]Gray D,Goodway B and Chen T.Bridging the gap：Using AVO to detect changes in fundamental elastic constants.SEG Technical Program Expanded Abstracts,1999,18：852-855.

[2]Downton J E,Lines L R.Constrained three parameter AVO inversion and uncertainty analysis.SEG Technical Program Expanded Abstracts,2001,20：251-254.

[3]Russell B H,Hedlin K,Hilterman F J et al.Fluid-property discrimination with AVO：A Biot-Gassmann perspective.Geophysics,2003,68(1)：29-39.

[4]楊培杰,印興耀,張欣.疊后地震盲反演及其應(yīng)用.石油地球物理勘探,2008,43(3)：284-290.

Yang Peijie,Yin Xingyao and Zhang Xin.Poststack seismic blind inversion and application.OGP,2008,43(3)：284-290.

[5]趙小龍,吳國忱,曹丹平.多尺度地震資料稀疏貝葉斯聯(lián)合反演方法.石油地球物理勘探,2016,51(6)：1156-1163.

Zhao Xiaolong,Wu Guochen,Cao Danping.A sparse Bayesian joint inversion of multi-scale seismic data.OGP,2016,51(6)：1156-1163.

[6]張豐麒,金之鈞,盛秀杰等.一種穩(wěn)健的彈性阻抗反演方法.石油地球物理勘探,2017,52(2)：294-303.

Zhang Fengqi,Jin Zhijun,Sheng Xiujie et al.A stable elastic impedance inversion approach.OGP,2017,52(2)：294-303.

[7]朱勇,唐斌.一種提高波阻抗反演建模精度的方法.江漢石油職工大學(xué)學(xué)報(bào),2017,20(6)：22-26.

Zhu Yong and Tang Bin.A method of improving modeling precision of wave impedance inversion.Journal of Jianghan Petroleum University of Staff and Workers,2017,20(6)：22-26.

[8]夏洪瑞,周開明,黃橋.波阻抗反演中的一種建模方法.石油物探,2004,43(1)：30-32.

Xia Hongrui,Zhou Kaiming and Huang Qiao.Building of initial model for wave impedance inversion.GPP,2004,43(1)：30-32.

[9]王偉,郭宇宏,高靜懷等.基于高精度地層傾角線性反演的地震全體積拉平方法.地球物理學(xué)報(bào),2013,56(3)：1012-1022.

Wang Wei,Guo Yuhong,Gao Jinghuai et al.Volumetric flattening through linear inversion scheme with accurate seismic dips.Chinese Journal of Geophysics,2013,56(3)：1012-1022.

[10]Schultz P S,Ronen S,Hattori M et al.Seismic-guided estimation of log properties,part 1：A data-driven interpretation methodology.The Leading Edge,1994,13(5)：305-310.

[11]Schultz P S,Ronen S,Hattori M et al.Seismic-guided estimation of log properties,part 3：A controlled study.The Leading Edge,1994,13(7)：770-776.

[12]Ronen S,Schultz P S,Hattori M et al.Seismic-guided estimation of log properties,part 2： Using artificial neural networks for nonlinear attribute calibration.The Leading Edge,1994,13(6)：674-678.

[13]Hampson D P,Schuelke J S and Quirein J A.Use of multi-attribute transforms to predict log properties from seismic data.Geophysics,2001,66(1)：220-236.

[14]Hale D.Image-guided blended neighbor interpolation of scattered data.SEG Technical Program Expanded Abstracts, 2009,28：1127-1130.

[15]Hale D.Image-guided 3D interpolation of borehole data.SEG Technical Program Expanded Abstracts,2010,29：1266-1270.

[16]劉賢紅,査樹貴.地震屬性約束初始波阻抗模型建立方法.中國石油學(xué)會(huì)2010年物探技術(shù)研討會(huì)論文集,2010,374-376.

Liu Xianhong and Zha Shugui.Seismic attribute constrained initial wave impedance modeling method.Symposium on Geophysical Prospecting Techniques,Chinese Petroleum Society,2010,374-376.

[17]Fomel S.Predictive painting of 3D seismic volumes.Geophysics,2010,75(4)：A25-A30.

[18]Karimi P and Fomel S.Image guided well log interpolation using predictive painting.SEG Technical Program Expanded Abstracts,2015,34：2786-2790.

[19]Fomel S.Applications of plane-wave destruction filters：Geophysics,2002,67(6)：1946-1960.

[20]Chen Z,Fomel S and Lu W.Accelerated plane-wave destruction.Geophysics,2013,78(1)：V1-V9.

[21]Claerbout J F.Earth Soundings Analysis：Processing Versus Inversion.Blackwell Scientific Publications Inc,1992.

[22]Fomel S.Shaping regularization in geophysical estimation problems.Geophysics,2007,72(2)：R29-R36.

[23]黃建平,李振春,孔雪等.基于PWD的繞射波波場分離成像方法綜述.地球物理學(xué)進(jìn)展,2012,27(6)：2499-2510.

Huang Jianping,Li Zhenchun,Kong Xue et al.The review of the wave field separation method about reflection and diffraction based on the PWD.Progress in Geophysics,2012,27(6)：2499-2510.

[24]Gan S,Wang S,Chen Y et al.Deblending using a structural-oriented median filter.SEG Technical Program Expanded Abstracts,2015,34：59-64.

[25]Tikhonov A N and Arsenin V Y.Solutions of Ill-posed Problems.New York：John-Wiley & Sons, 1977.

[26]王治強(qiáng),曹思遠(yuǎn),陳紅靈等.基于TV約束和Toepplitz矩陣分解的波阻抗反演.石油地球物理勘探,2017,52(6)：1193-1199.

Wang Zhiqiang,Cao Siyuan,Chen Hongling et al.Wave impedance inversion based on TV regularization and Toeplitz-sparse matrix factorization.OGP,2017,52(6)：1193-1199.