王治強 曹思遠 陳紅靈 孫曉明 樊 平

(①中國石油大學(xué)地球物理與信息工程學(xué)院,北京102249； ②中國石油大學(xué)油氣資源與探測國家重點實驗室,北京102249; ③東方地球物理公司遼河物探處,遼寧盤錦124000)

·偏移成像·

基于TV約束和Toeplitz矩陣分解的波阻抗反演

王治強*①②曹思遠①②陳紅靈①②孫曉明①②樊 平③

(①中國石油大學(xué)地球物理與信息工程學(xué)院,北京102249； ②中國石油大學(xué)油氣資源與探測國家重點實驗室,北京102249; ③東方地球物理公司遼河物探處,遼寧盤錦124000)

利用波阻抗剖面的非高斯分布特點以及地震子波褶積矩陣的Toeplitz結(jié)構(gòu)，對波阻抗剖面進行全變分(TV)約束，可以在壓制隨機噪聲的同時保持剖面的不連續(xù)性，對地震子波褶積矩陣進行Toeplitz稀疏矩陣分解得到地震子波的稀疏表達。地震資料的低頻損失導(dǎo)致無法反演出波阻抗的低頻背景，故將測井或解釋層位信息通過最小二乘法建立約束條件。建立的約束最優(yōu)化目標(biāo)函數(shù)可以同時反演子波和波阻抗。文中將地震剖面整體處理，反演結(jié)果比常規(guī)的逐道反演具有更高的精度、橫向連續(xù)性和抗噪性。

波阻抗 TV約束 Toeplitz結(jié)構(gòu) 矩陣分解 約束最優(yōu)化

1 引言

波阻抗反演是疊后地震資料解釋中非常重要的技術(shù)之一。它將地震剖面轉(zhuǎn)換為波阻抗剖面，不僅便于解釋人員直接對比地震資料與測井資料，而且能對地層物性參數(shù)的變化進行有效分析，得到其在空間的分布規(guī)律，指導(dǎo)油氣勘探開發(fā)[1]。常規(guī)波阻抗反演方法是在反褶積的基礎(chǔ)上對提取的反射系數(shù)逐道進行遞推反演[2]。地震資料的低頻損失導(dǎo)致無法反演出波阻抗的低頻背景[3]，故將測井或解釋層位信息通過最小二乘法建立約束條件。隨著石油勘探開發(fā)的不斷深入，基于同相軸追蹤進行層位解釋的常規(guī)方法已經(jīng)無法滿足目前的地質(zhì)需求。借助于測井資料和地質(zhì)認識，利用波阻抗反演資料綜合地震資料識別地下構(gòu)造、巖性已經(jīng)成為地震資料解釋和儲層預(yù)測的重要技術(shù)[4]。然而，地震資料低頻和高頻信息的缺失、地震子波的未知、噪聲的干擾以及地質(zhì)結(jié)構(gòu)的復(fù)雜性降低了地震反演的精度。根據(jù)正演模型的不同，波阻抗反演分為基于波動方程的反演和基于褶積模型的反演兩大類。前者算法結(jié)構(gòu)復(fù)雜、計算量大、抗干擾性差，很難獲得一個穩(wěn)定解，未得到廣泛應(yīng)用。目前主要使用基于褶積模型的反演方法，該算法相對簡單，對噪聲敏感性相對較小[5]。

波阻抗反演技術(shù)經(jīng)歷了從直接遞推反演到基于模型的迭代反演的發(fā)展過程[6]。直接遞推反演分兩步完成，首先對地震資料進行反褶積處理，求得反射系數(shù)剖面，再對反射系數(shù)序列進行道積分得到波阻抗剖面。

基于模型的迭代反演方法以最小二乘法為基礎(chǔ)，約束正演模型和實際監(jiān)測數(shù)據(jù)之間的誤差，再根據(jù)地下波阻抗的物性特征加入約束條件，建立約束最優(yōu)化目標(biāo)函數(shù)。通過不斷修改初始模型使目標(biāo)函數(shù)取得極值，即可得到最合理的解。迭代反演能夠有效克服地震資料直接反演存在的帶限問題[7]。

根據(jù)模型與觀測數(shù)據(jù)之間的函數(shù)關(guān)系，反演問題分為線性反演和非線性反演。Cooke等[8]將波阻抗反演問題線性化，使用廣義線性反演方法直接從地震剖面中逐道估計波阻抗； Ma[9]在假定地震子波已知的前提下，使用模擬退火算法實現(xiàn)了單道非線性波阻抗反演； Hamid等[10]同樣在地震子波已知的假設(shè)條件下，在對數(shù)域?qū)Ψ囱輪栴}線性化，使用Tikhonov正則化方法同時對多道數(shù)據(jù)進行波阻抗反演。實際上，Tikhonov正則化方法使目標(biāo)函數(shù)的解過度光滑，無法滿足地質(zhì)構(gòu)造復(fù)雜的情況[11]。

本文充分利用地下巖石波阻抗的分布特征、地震子波褶積矩陣的Toeplitz結(jié)構(gòu)以及地震子波的光滑性建立約束條件。全變分(TV)約束可以在壓制隨機噪聲的同時很好地保持地層的邊界特征[3]。對子波褶積矩陣進行稀疏矩陣分解可以得到地震子波的稀疏表達[12]。由于地震資料的低頻損失，無法反演出波阻抗的低頻背景，故將測井或解釋層位信息通過最小二乘法建立約束條件[3]。建立的約束最優(yōu)化目標(biāo)函數(shù)通過基于褶積模型的迭代反演方法求解波阻抗。

由于建立的約束最優(yōu)化目標(biāo)函數(shù)是非凸的，求解困難，文中采用了對阻抗剖面和子波交替更新求解的策略，將非凸優(yōu)化問題退化為兩個凸優(yōu)化問題分別求解，有效降低了求解難度[13]。本文將地震剖面整體處理，反演結(jié)果比常規(guī)的逐道反演具有更好的橫向連續(xù)性和抗噪能力。

2 波阻抗的全變分(TV)正則化約束

在二維條件下，假設(shè)反射系數(shù)剖面和波阻抗剖面分別為

(1)

R、Z滿足

(2)

可引入變量

(3)

則有

R=L1X

(4)

其中

(5)

因此，地震記錄的褶積模型可表示為

D=WL1X+N

(6)

式中：D為地震記錄；W為n×n維地震子波褶積矩陣;N為n×m維加性噪聲。

波阻抗反演的目的就是在地震記錄D已知的條件下求解X，再通過指數(shù)運算得到波阻抗剖面Z。最直接的方法是直接求解式(6)中的X

即

(7)

式(7)是一個病態(tài)問題，無法得到適定的解[14]。為得到適定解，對X進行全變分正則化約束，全變分約束可以在壓制噪聲的同時保持阻抗剖面邊界的不連續(xù)性[3]，得到

(8)

式中：ε為噪聲水平； ‖D-WL1X‖2F≤ε是最基本的最小二乘約束； F為范數(shù)； ‖X‖TV是全變分約束。

地震資料的低頻損失導(dǎo)致無法反演出波阻抗的低頻背景，故將測井或解釋層位信息以最小二乘約束‖L2X-Xlow‖2F作為約束條件，L2為低通濾波算子，Xlow為常用對數(shù)阻抗低頻背景。引入拉格朗日乘子λ1、λ2，將約束最優(yōu)化問題轉(zhuǎn)換為無約束最優(yōu)化問題，建立目標(biāo)函數(shù)

‖D-WL1X‖2F+λ1‖X‖TV+

λ2‖L2X-Xlow‖2F

(9)

上述波阻抗反演方法的不足之處在于首先需要提取比較精確的子波，步驟繁瑣[15]。

3 Toeplitz結(jié)構(gòu)的稀疏矩陣分解

子波褶積矩陣W具有Toeplitz矩陣結(jié)構(gòu)[12]，對于一般的W做如下分解

即

W=wn-1In-1+…+w1I1+w0I0+w-1I-1+

(10)

式中矩陣I為一系列相應(yīng)的斜對角矩陣。令向量a=[wn-1,wn-2,…,w0,…,w-(n-2),w-(n-1)]T

目標(biāo)函數(shù)可修改為

λ1‖X‖TV+λ2‖L2X-Xlow‖2F

(11)

實際地震子波長度l遠遠小于記錄長度n，即(l?n)。地震子波的褶積矩陣為[12]

(12)

對于向量a而言，只有w0,w1,…,wl-1為非零值，其他元素均為0。表明向量a同樣具有稀疏性，可以用L1范數(shù)作為約束條件。約束最優(yōu)化目標(biāo)函數(shù)改進為

λ2‖L2X-Xlow‖2F+λ3‖a‖1

(13)

考慮到地震子波在一定程度上還具有光滑性[12]，依然可以作為反演的約束條件。光滑性可以用差分系數(shù)的稀疏性來衡量，給目標(biāo)函數(shù)加入以下形式的懲罰項

|ak-ak-1|=λ4‖L3a‖1

(14)

其中

L3[(2n-2)×(2n-1)]

目標(biāo)函數(shù)可改進為

λ2‖L2X-Xlow‖2F+λ3‖a‖1+λ4‖L3a‖1

(15)

求解式(15)是極其困難的，因為該問題屬于非凸優(yōu)化問題。如果固定a或者X中的其中一個，則非凸優(yōu)化問題退化為一個凸優(yōu)化問題?？梢圆捎媒惶婀潭ㄗ兞康姆椒ㄖ鸩降蠼庠瓎栴}[12]。

首先，固定子波褶積矩陣W，問題退化為

‖D-WL1X‖2F+λ1‖X‖TV+

λ2‖L2X-Xlow‖2F

(16)

同樣，若固定X，則問題退化為

λ4‖L2a‖1

(17)

為了方便求解式(17)，將式(17)整理為

(18)

綜上所述，通過不斷交替更新子波矩陣W和X，直到得到滿意的結(jié)果。最開始固定的是W，需要一個初始子波。本文方法在大多數(shù)情況下直接以地震資料頻譜的質(zhì)心頻率為主頻的零相位雷克子波作為初始子波就可以得到滿意的結(jié)果。

上述過程得到對數(shù)波阻抗剖面X，利用式(19)可得到波阻抗剖面

Z(i,j)=e2X(i,j)

(19)

4 模型測試

圖1 不含隨機噪聲模型反演結(jié)果

為了測試反演方法的抗噪性，在合成記錄中加入5%的隨機噪聲得到正演地震剖面。圖2為含隨機噪聲5%模型反演結(jié)果?？梢钥闯?，與無噪聲反演結(jié)果對比，阻抗剖面的橫向連續(xù)性較強，但阻抗值在層內(nèi)出現(xiàn)微小的擾動(尤其圖2c中的第5層和第7層)。但對于阻抗剖面整體而言，這種微小擾動完全不影響判斷地下巖性及其空間展布。反演地震子波與實際地震子波波形差別較小(圖2e)。

圖2 含隨機噪聲5%模型反演結(jié)果

圖3 含隨機噪聲10%模型反演結(jié)果

上述結(jié)果表明，在弱噪聲條件下，文中的反演方法可以得到相對精確的阻抗剖面。合成記錄加入10%的隨機噪聲得到地震剖面。圖3為含隨機噪聲10%模型反演結(jié)果?？梢钥闯觯瑢?shù)阻抗剖面依然保持相對良好的橫向連續(xù)性，但阻抗值受噪聲影響較大，反演結(jié)果與給定模型(圖1a)有顯著區(qū)別，但阻抗之間的相對大小并未發(fā)生太大變化。由圖3d可見，反演的子波與實際波形有較大差別，這是由于受隨機噪聲的影響。

綜上所述，通過反演含有不同強度隨機噪聲的模型發(fā)現(xiàn)，本文方法在隨機噪聲強度小于約6%的情況下，可以反演出相對精確的阻抗剖面，提取的地震子波波形與給定子波一致性較好。當(dāng)隨機噪聲強度接近10%時，反演阻抗剖面的局部穩(wěn)定性變差，但反演結(jié)果保持了較好的橫向連續(xù)性，依然可以提取到相對準確的地震子波，反演阻抗剖面可以反映地下巖性及其空間展布。

5 實際資料處理

圖4是A區(qū)實際資料疊加剖面，該剖面過KE031井和JI14井，假設(shè)地震資料已經(jīng)實現(xiàn)吸收衰減等補償(Q=∞)[18],分別截取黃色方框內(nèi)剖面作為本次反演的原始數(shù)據(jù)。圖5是圖4中黃色區(qū)域內(nèi)截取的地震剖面(圖5a、圖5b分別是KE031井和JI14井的過井剖面)，分別提取測井曲線的低頻部分作為反演的低頻背景[19](圖6)。圖7為圖5數(shù)據(jù)不同方法反演的波阻抗剖面。可以看出，本文方法反演的波阻抗剖面界面連續(xù)性更好，薄層分辨率更高(圖7c、圖7d)。為了驗證方法的可靠性，分別將KE031井和JI14井的測井曲線標(biāo)定在反演剖面上，如圖7中的黑色曲線所示。可以看出，本文方法反演的波阻抗剖面與測井曲線更加吻合，可以提高地下巖性和薄層的識辨能力[20-22]。圖8為本文方法同步反演的地震子波，初始子波為主頻為28Hz(根據(jù)地震資料確定)的零相位雷克子波?？梢钥闯觯囱葑硬ㄅc初始子波波形基本吻合。

圖6 圖5數(shù)據(jù)的低頻背景

圖7 圖5a(左)及圖5b(右)數(shù)據(jù)不同方法反演結(jié)果對比

圖8 本文方法反演的地震子波

6 結(jié)論

本文提出的基于TV約束和Toeplitz矩陣分解的波阻抗反演與常規(guī)波阻抗反演相比，主要有以下幾點創(chuàng)新和優(yōu)勢。

(1)將子波提取與波阻抗反演合并為一個問題，建立并求解約束最優(yōu)化目標(biāo)函數(shù)，直接從地震資料中求取地震子波和波阻抗。

(2)利用地震子波褶積矩陣的Toeplitz結(jié)構(gòu)特征，運用稀疏矩陣分解[12]對原地震子波的褶積矩陣進行分解及重構(gòu)，將地震子波的光滑性及連續(xù)性作為約束條件。

(3)利用全變分約束(TV)可以在壓制隨機噪聲的同時很好刻畫阻抗剖面的邊界特性，保持其不連續(xù)性的特點。由于地震資料的低頻損失，直接反演出的結(jié)果只代表實際阻抗的中頻成分，文中設(shè)計濾波算子利用最小二乘原理將測井信息作為約束條件。

(4)直接求解建立的目標(biāo)函數(shù)為非凸優(yōu)化問題非常困難。通過交替固定地震子波和對數(shù)阻抗剖面，將一個非凸優(yōu)化問題退化為兩個凸優(yōu)化子問題，降低了求解難度。

本文方法的不足之處在于該反演方法基于多道處理，受計算機內(nèi)存限制，無法同時進行較大數(shù)據(jù)量的運算。

[1] 撒利明,楊午陽,姚逢昌等.地震反演技術(shù)回顧與展望.石油地球物理勘探,2015,50(1):184-202.

Sa Liming,Yang Wuyang,Yao Fengchang et al.Review and expectation of seismic inversion technology.OGP,2015,50(1):184-202.

[2] 呂鐵良.波阻抗約束反演中的約束方法研究[學(xué)位論文].山東青島：中國石油大學(xué),2007.

Lü Tieliang.Constraint Method Research in Impe-dance Inversion[D].China University of Petroleum,Qingdao,Shandong，2007.

[3] Gholami A.Nonlinear multichannel impedance inversion by total-variation regularization.Geophysics,2015,80(5)：R217-R224.

[4] 彭傳平.寬帶約束反演方法實現(xiàn)及應(yīng)用研究[學(xué)位論文].陜西西安：長安大學(xué),2008.

Peng Chuanzhen.Broadband Constraint Inversion Method and Application Research[D].Chang’an university,Xi’an,Shaanxi,2008.

[5] 王圣川.正則化約束稀疏脈沖地震反演方法及應(yīng)用研究[學(xué)位論文].四川成都：電子科技大學(xué),2014.

Wang Shengchuang.Regularization Constraint Sparse Pulse Seismic Inversion Method and Application Research[D].University of Electronic Science and Technology of China,Chengdu,Sichuan，2014.

[6] 王東燕.地震約束波阻抗反演及應(yīng)用研究[學(xué)位論文].陜西西安：長安大學(xué),2006.

Wang Dongyan.Impedance Inversion Based on Seismic Constraint and Application Research[D].Chang’an University,Xi’an,Shaanxi，2006.

[7] 裴云龍.稀疏約束反褶積及其波阻抗反演方法研究[學(xué)位論文].北京：中國地質(zhì)大學(xué)(北京),2009.

Pei Yunlong.Sparse Constraint Deconvolution and Wave Impedance Inversion Method Research[D].China University of Geosciences(Beijing),Beijing,2009.

[8] Cooke D A and Schneider W A.Generalized linear inversion of reflection seismic data.Geophysics,1983,48(6):665-676.

[9] Ma X Q.A constrained global inversion methon using an overparameterized scheme:Application to post-stack seismic data.Geophysics,2012,66(2):613-626.

[10] Hamid H and Pidlisecky A.Multitrace impedance inversion with lateral constraints.Geophysics,2015,80(6):M101-M111.

[11] Gholami A.A fast automatic multichannel blind seismic inversion for high-resolution impedance recovery.Geophysics,2016,81(5):V357-V364.

[12] Wang Lingling,Zhao Qian,Gao Jinghuai.Seismic sparse-spike deconvolution via Toeplotz-sparse matrix factorization.Geophysics,2016,81(3):V169-V182.

[13] 印興耀，曹丹平,王寶麗等.基于疊前地震反演的流體識別方法研究進展.石油地球物理勘探,2014,49(1):22-34.

Yin Xingyao,Cao Danping,Wang Baoli et al.Research progress of fluid discrimination with prestack seismic inversion.OGP,2014,49(1):22-34.

[14] Gholami A,Sacchi M D.Fast 3D blind seismic deconvolution via constrained total variation and GCV.Society for Industrial and Applied Mathema-tics,2013,6(4):2350-2369.

[15] 孫學(xué)凱,馮世民.地震反演系統(tǒng)中的子波提取方法.物探化探計算技術(shù),2010,32(2):120-125.

Sun Xuekai,Feng Shimin.Wavelet extraction methods of seismic inversion system.Geophysical Computing Technology,2010,32(2):120-125.

[16] 朱衛(wèi)星,楊玉卿,趙永生等.測井地震聯(lián)合反演在地質(zhì)導(dǎo)向風(fēng)險控制中的應(yīng)用.石油地球物理勘探,2013,48(增刊1):181-185.

Zhu Weixing,Yang Yuqing,Zhao Yongsheng et al.Joint inversion of logging and seismic data for geosteering risk control.OGP,2013,48(S1):181-185.

[17] Dai Ronghuo,Zhang Fanchang and Liu Hanqing.Seismic inversion based on proximal objective function optimization algorithm.Geophysics,2016,81(5):R237-R246.

[18] 吳華.地震儲層反演方法研究及應(yīng)用[學(xué)位論文].北京：中國地質(zhì)大學(xué)(北京),2015.

Wu Hua.Research and Application of Seismic Reservoir Inversion method[D].China University of Geosciences(Beijing),Beijing,2015.

[19] 欒穎.約束稀疏脈沖波阻抗反演方法在煤層識別中的應(yīng)用[學(xué)位論文].吉林長春：吉林大學(xué),2010.

Luan Ying.Constrained Sparse Pulse Wave Impe-dance Inversion Method in the Application of the Coal Seam Recognition[D].Jilin University,Changchun，Jilin,2010.

[20] 張國偉.疊后波阻抗反演及儲層預(yù)測[學(xué)位論文].四川成都：成都理工大學(xué),2014.

Zhang Guowei.Poststack Wave Impedance Inversion and Reservoir Prediction[D].Chengdu Univerisity of Technology,Chengdu，Sichuan，2014.

[21] 楊立強.測井約束地震反演綜述.地球物理學(xué)進展,2003,18(3):530-534.

Yang Liqiang.Well logging constrained seismic inversion reviewed.Progress in Geophysics,2003,18(3)：530-534.

[22] 石艷玲，黃文輝，魏強等.電磁井震約束反演識別川中深層裂谷.石油地球物理勘探,2016,51(6):1233-1240.

Shi Yanling,Huang Wenhui,Wei Qiang et al.Identify deep rift in central Sichuan via inversion with constraint of electromagnetic and well-seismic.OGP,2016,51(6):1233-1240.

*北京市昌平區(qū)府學(xué)路18號中國石油大學(xué)(北京)地球物理與信息工程學(xué)院，102249。Email:529433961@qq.com

本文于2016年10月28日收到，最終修改稿于2017年8月27日收到。

1000-7210(2017)06-1193-07

王治強，曹思遠，陳紅靈，孫曉明，樊平.基于TV約束和Toeplitz矩陣分解的波阻抗反演.石油地球物理勘探，2017,52(6):1193-1199,1245.

P631

A

10.13810/j.cnki.issn.1000-7210.2017.06.009

(本文編輯：金文昱)

王治強 碩士研究生，1991年生；2015年本科畢業(yè)于長安大學(xué)勘查技術(shù)與工程專業(yè)； 2015年至今在中國石油大學(xué)(北京)地球物理與信息工程學(xué)院攻讀地質(zhì)資源與地質(zhì)工程碩士學(xué)位。主要研究方向：地震資料處理及其疊后反演。