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

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

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

·偏移成像·

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

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

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

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

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

1 引言

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

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

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

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

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

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

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

在二維條件下，假設反射系數(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是全變分約束。

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

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

λ2‖L2X-Xlow‖2F

(9)

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

3 Toeplitz結構的稀疏矩陣分解

子波褶積矩陣W具有Toeplitz矩陣結構[12]，對于一般的W做如下分解

即

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

(10)

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

目標函數(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)化目標函數(shù)改進為

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

(13)

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

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

(14)

其中

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

目標函數(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，直到得到滿意的結果。最開始固定的是W，需要一個初始子波。本文方法在大多數(shù)情況下直接以地震資料頻譜的質心頻率為主頻的零相位雷克子波作為初始子波就可以得到滿意的結果。

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

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

(19)

4 模型測試

圖1 不含隨機噪聲模型反演結果

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

圖2 含隨機噪聲5%模型反演結果

圖3 含隨機噪聲10%模型反演結果

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

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

5 實際資料處理

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

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

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

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

6 結論

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

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

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

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

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

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

[1] 撒利明,楊午陽,姚逢昌等.地震反演技術回顧與展望.石油地球物理勘探,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] 呂鐵良.波阻抗約束反演中的約束方法研究[學位論文].山東青島：中國石油大學,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)及應用研究[學位論文].陜西西安：長安大學,2008.

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

[5] 王圣川.正則化約束稀疏脈沖地震反演方法及應用研究[學位論文].四川成都：電子科技大學,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] 王東燕.地震約束波阻抗反演及應用研究[學位論文].陜西西安：長安大學,2006.

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

[7] 裴云龍.稀疏約束反褶積及其波阻抗反演方法研究[學位論文].北京：中國地質大學(北京),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] 孫學凱,馮世民.地震反演系統(tǒng)中的子波提取方法.物探化探計算技術,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)合反演在地質導向風險控制中的應用.石油地球物理勘探,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] 吳華.地震儲層反演方法研究及應用[學位論文].北京：中國地質大學(北京),2015.

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

[19] 欒穎.約束稀疏脈沖波阻抗反演方法在煤層識別中的應用[學位論文].吉林長春：吉林大學,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] 張國偉.疊后波阻抗反演及儲層預測[學位論文].四川成都：成都理工大學,2014.

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

[21] 楊立強.測井約束地震反演綜述.地球物理學進展,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ū)府學路18號中國石油大學(北京)地球物理與信息工程學院，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è)于長安大學勘查技術與工程專業(yè)； 2015年至今在中國石油大學(北京)地球物理與信息工程學院攻讀地質資源與地質工程碩士學位。主要研究方向：地震資料處理及其疊后反演。