基于平面波解構的初始波阻抗建模方法

李海山　楊午陽　雍學善

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

1　引言

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

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

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

2　基于PWD的局部斜率計算原理

2.1　平面波解構算子

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

(1)

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

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

(2)

其中f(t)代表波場。若斜率σ不依賴于時間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ù)可以通過在零頻率附近進行泰勒級數(shù)展開確定，其中通常采用的三點中心濾波器B3(zt)為[19]

(8)

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

(9)

2.2　局部斜率計算方法

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

C(σ)u≈0

(10)

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

(11)

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

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

(12)

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

εDΔσ≈0

(13)

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

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

3.1　初始波阻抗建模原理

計算得到同相軸局部斜率之后，就可以利用平面波解構濾波器通過已知道的波阻抗值預測相鄰道的波阻抗值。如果采用式(9)給出的三點中心濾波器系數(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ù)(1≤t≤T)，系數(shù)

(19)

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

(20)

如果由x道預測x+1道，則

(21)

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

3.2　Tikhonov正則化

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

Cu≈w

(22)

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

(23)

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

(24)

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

(25)

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

(26)

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

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

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

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

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

5　實際資料應用

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

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

圖2　實際資料應用結果比較

6　結論

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

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

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

[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]楊培杰,印興耀,張欣.疊后地震盲反演及其應用.石油地球物理勘探,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]朱勇,唐斌.一種提高波阻抗反演建模精度的方法.江漢石油職工大學學報,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]王偉,郭宇宏,高靜懷等.基于高精度地層傾角線性反演的地震全體積拉平方法.地球物理學報,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]劉賢紅,査樹貴.地震屬性約束初始波阻抗模型建立方法.中國石油學會2010年物探技術研討會論文集,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的繞射波波場分離成像方法綜述.地球物理學進展,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]王治強,曹思遠,陳紅靈等.基于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.