非均質HTI介質裂縫弱度參數地震散射反演

潘新朋 張廣智 印興耀

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580；②海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266071)

·綜合研究·

非均質HTI介質裂縫弱度參數地震散射反演

潘新朋*①張廣智①②印興耀①②

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580；②海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266071)

長波長假設下各向同性地層中發育的垂直排列裂縫可構成等效的具有水平對稱軸的橫向各向同性(HTI)介質。裂縫弱度是重要的非均質各向異性介質特征參數，彈性逆散射理論是非均質介質參數反演的有效途徑?；诘卣鹕⑸淅碚?，首先推導了非均質HTI介質中含裂縫弱度參數的縱波散射系數方程，并通過裂縫弱度的參數組合，提出了一種新穎的HTI介質方位彈性阻抗(AEI)參數化方法。為了提高反演的穩定性和橫向連續性，改進了貝葉斯框架下的HTI介質各向異性方位彈性阻抗反演方法，同時考慮了柯西稀疏約束正則化和平滑模型約束正則化，最終運用非線性迭代重加權最小二乘策略實現了裂縫弱度參數的穩定估算。模型和實際資料處理表明，反演結果與測井解釋數據相吻合，證明利用該方法能穩定可靠地從方位地震資料中反演得到裂縫弱度參數。

HTI各向異性 裂縫弱度 彈性逆散射 方位彈性 阻抗正則化

1 引言

非均質HTI介質中地震波速度與波傳播方向相關，即存在地震各向異性現象[1，2]。在長波長假設條件下，各向同性背景地層中發育一組垂直排列的裂縫可等效為具有水平對稱軸的橫向各向同性(HTI)介質[3，4]。連通的裂縫不僅能為油氣儲存提供孔隙空間，而且增加了巖石滲透率。因此，利用巖石物理分析、地質、測井及地震資料綜合預測地下裂縫，識別裂縫型儲層的位置，估算裂縫特征參數，實現裂縫特征的有效描述，對裂縫儲層的勘探和油氣開發具有現實意義[5，6]。裂縫弱度參數是裂縫特征描述的重要參數，探索裂縫弱度參數與地震反射之間的關系尤為關鍵[7，8]。

由于非均質HTI介質中地震波傳播的復雜性，很難推導其線性的縱波反射系數解析表達式。基于弱各向異性理論[9]，Tsvankin[10]引入了表征HTI介質的縱波各向異性符號。在此基礎上，前人推導了HTI介質中的縱波反射系數公式[3，4，11，12]。然而，這些推導過程受限于水平層狀假設，當地下存在裂縫等非均質介質時，就會產生散射[13]，此時水平層狀假設不再滿足非均質各向異性介質參數的反演需求。彈性逆散射理論是非均質介質參數反演的有效途徑[14，15]，如常用的Born或Rytov近似[16，17]?；趶椥阅嫔⑸淅碚?，Shaw等[18，19]借助穩相法推導了非均質各向異性介質的縱波散射系數方程。據此，本文也應用彈性逆散射原理推導了由裂縫弱度參數表征的HTI介質縱波散射系數方程。

由于分方位地震數據的低信噪比，各向異性參數的反演仍然是一個難題。由于角度疊加數據具較高信噪比，彈性阻抗(EI)反演已經廣泛應用于各向同性介質中[20-24]，目前也逐步擴展到各向異性介質[25-28]。然而，所推導并應用的EI方程的各向異性擾動部分存在指數校正項，與各向同性背景部分的表達式存在明顯差異，可能增加各向異性參數反演的不穩定性，導致各向異性參數反演的精度低于各向同性參數反演。為此，本文提出一種新的HTI介質各向異性參數化方法，并推導了相應的各向異性方位彈性阻抗方程。為了提高反演的穩定性和橫向連續性，發展了貝葉斯框架下的各向異性方位彈性阻抗反演方法，同時考慮了柯西稀疏約束正則化和平滑模型約束正則化[29-32]，應用非線性迭代重加權最小二乘策略[33]實現了各向異性參數的穩定估算。模型和實際資料處理表明，反演結果與測井解釋數據相吻合，證明該方法能穩定可靠地從方位疊前地震資料中獲取裂縫弱度參數，為非均質HTI介質各向異性預測提供了一種可靠穩定的地震反演途徑。

2 方法原理

2.1 非均質HTI介質有效彈性剛度矩陣

基于Schoenberg線性滑移理論[34，35]，HTI介質有效剛度矩陣SHTI可簡單地由各向同性背景柔度矩陣Sb與附加的裂縫柔度矩陣Sf之和表征，即

SHTI=Sb+Sf

(1)

因此，各向同性背景中發育一組平行的垂直裂縫的等效HTI介質的有效彈性剛度矩陣CHTI可表示為

(2)

裂縫弱度參數是儲層裂縫描述的特征參數，有助于指導地下的裂縫識別。 Hsu等[7]引入了一組無量綱的描述裂縫的特征參數

(3)

式中： 下標N、V和H分別表示法向、垂直切向和水平切向； 裂縫柔度矩陣元素KN、KV、KH、KNV、KNH和KVH均表征線性滑移邊界條件中的附加裂縫柔度；Mb和μb是背景介質的縱波模量和橫波模量；ΔN、ΔV和ΔH分別表征法向、垂直切向及水平切向裂縫弱度； 其他參數類似。因此，HTI介質的有效彈性剛度矩陣CHTI可用裂縫弱度參數表征為

(4)

其中

(5)

(6)

根據一階擾動理論，HTI介質有效剛度矩陣可表示為均質各向同性背景矩陣與其擾動矩陣之和

(7)

式中： ΔM和Δμ分別表示界面兩側縱波模量與橫波模量的差值； δΔN、δΔV和δΔH等參數分別表示界面兩側裂縫弱度參數的差值。

2.2 非均質各向異性HTI介質縱波散射系數

在弱散射假設條件下，基于Born近似和穩相法，非均質HTI介質縱波散射系數可表示為[18，19，36]

(8)

式中：θ是入射角；r表示激發點到接收點距離；r0是水平界面一點；S(r0)是散射函數，可表示為

S(r0)=Δρξ+ΔCηmn

(9)

其中

(10)

式中p和t分別表示慢度矢量和極化矢量。

在此基礎上，推導了非均質HTI介質的縱波散射系數方程(詳見附錄A)

=a(θ)RM+b(θ)Rμ+c(θ)Rρ+

(11)

(12)

實際儲層裂縫大多具有旋轉不變特征屬性，其剛度張量矩陣滿足以下條件

(13)

(14)

其中

h(θ,φ)=2g(sin2θcos2φ-sin2θtan2θsin2φcos2φ)

(15)

(16)

(17)

式中

2.3 非均質HTI介質方位彈性阻抗參數化

Connolly[20]定義彈性阻抗為

(18)

式中：n-1和n分別代表下層和上層介質； EI表示介質的彈性阻抗屬性。對于方位各向異性介質而言，RPP≡RPP(θ,φ)且EI≡EI(θ,φ)。

單界面的縱波散射系數方程(式(18))可擴展為時間連續的散射系數函數[37]

(19)

對式(19)進行積分計算，可得

lnEI(t,θ,φ)=a(t,θ)lnM(t)+b(t,θ)lnμ(t)+

(20)

以指數形式表示式 (20)，則非均質HTI介質的彈性阻抗可表示為

EI(t,θ,φ)=[M(t)]a(t,θ)·[μ(t)]b(t,θ)·[ρ(t)]c(t,θ)·

(21)

為了消除式(21)中入射角變化對量綱尺度的影響，對其進行歸一化[21]，可得

(22)

2.4 非均質HTI介質方位彈性阻抗(AEI)反演

非均質HTI介質的AEI反演主要包括方位各向異性彈性阻抗數據體的反演及介質彈性參數和各向異性特征參數的提取。不同方位、不同角度的各向異性彈性阻抗數據體可使用常用的疊后反演方法估算，如常用的稀疏約束脈沖反演(CSSI)方法，但各向異性參數的反演到目前為止仍是一個難題。為了提高參數反演的穩定性和抗噪性，我們發展了一種貝葉斯框架下的非均質HTI介質AEI穩定反演方法，同時考慮柯西稀疏約束正則化和平滑模型約束正則化，最終使用非線性的迭代重加權最小二乘策略實現了各向異性參數的穩定估算。

對式(22)兩邊取對數，可得到用于非均質HTI介質彈性參數和各向異性參數估算的線性化AEI表達式

(23)

該式寫成矩陣形式為

LEI=Gm

(24)

其中

(25a)

(25b)

對式(24)中的參數進行去相關，得

LEI=G′m′

(26)

式中G′和m′分別表示去相關后的系數矩陣及對數域參數的散射系數。

假定地震噪聲服從高斯分布，而待反演參數服從修正柯西分布[29-32]，則后驗概率密度函數可表示為

(27)

地震數據的低頻分量相對較弱，易受噪聲影響。平滑模型約束正則化能夠補充反演參數的低頻信息，提高其穩定性及橫向連續性[31，32]。對式(27)兩邊取對數，并添加平滑模型約束正則化項，得

F(m′)=(LEI-G′m′)T(LEI-G′m′)+

(28)

其中

(29)

式中：λCauchy表示修正柯西稀疏約束正則化因子；λi為彈性參數與各向異性特征參數平滑模型約束的正則化參數；P表示積分算子矩陣。

柯西分布約束正則化項及平滑模型約束正則化項的引入使得方程式(28)具有弱非線性，采用重加權最小二乘反演方法(IRLS)能夠穩定的獲得參數的反演結果[33]。

3 模型試算與實際資料處理

3.1 模型試算

基于HTI介質巖石物理等效模型，可使用常規測井資料計算HTI介質的彈性參數及裂縫各向異性參數，從而建立參數的背景平滑模型，構建模型的具體步驟可參考文獻[38]。

基于HTI介質巖石物理等效模型，本文首先使用測井解釋數據估測了HTI介質的彈性參數及裂縫各向異性特征參數。為了驗證非均質HTI介質各向異性參數地震散射反演方法的可行性和穩定性，本文采用實際的單井模型進行驗證。模型測試使用主頻35Hz的Ricker子波。分別使用精確Zoeppritz方程與褶積模型正演合成的無噪及信噪比(SNR)分別為4和2的角度域疊前方位道集進行方法的驗證。圖1分別是無噪及含噪方位疊前道集反演獲得的參數結果，其中黑色曲線為測井真實數據，綠色曲線為初始背景模型，紅色曲線為反演結果。反演結果表明，在無噪聲情況下，該方法能夠獲取與真實值基本吻合的縱、橫波模量、密度及裂縫各向異性參數，驗證了本文方法的可行性。在適當噪聲的情況下，縱、橫波模量、密度及裂縫各向異性參數的反演結果與真實值仍有較高的吻合度。

3.2 實際資料處理

實際方位地震道集來自于中國陸上某勘探工區，進行疊前反演前，需對地震數據進行保幅、去噪等預處理。筆者利用本文提出的方法對不同角度疊加剖面進行了實際資料的試處理。本文反演方法輸入的是不同方位的部分角度疊加道集，圖2a、圖2b及圖2c分別是基于CSSI反演的AEI剖面； 圖3、圖4分別是估算獲得的縱波模量、橫波模量、密度彈性參數及裂縫各向異性參數剖面。同時，各圖中分別添加了對應的測井曲線，與地震反演結果進行比照。

從圖中可以看出，利用反演的AEI數據體提取的縱波模量、橫波模量及密度參數和裂縫各向異性參數結果與測井數值之間吻合較好，變化趨勢保持一致。同時，估測的裂縫各向異性參數值在裂縫發育位置呈現高值特征，與實際鉆遇結果吻合較好(紅色虛線橢圓為裂縫型油藏發育區)。因此，在實際反演得到裂縫各向異性參數后， 需要結合測井曲線及巖石物理分析結果，圈定裂縫各向異性參數值升高的區域為裂縫發育區域。

圖1 無噪及含噪情況下彈性參數及裂縫各向異性參數反演結果

圖2 不同方位的小、中、大角度彈性阻抗反演結果

圖3 彈性參數反演結果

圖4 裂縫各向異性參數反演結果

4 結論

裂縫弱度參數的可靠估算有助于非均質HTI介質的特征描述，利用AEI反演獲取裂縫弱度參數成為利用方位地震資料進行裂縫發育區識別的重要手段。地震散射理論是非均質介質參數反演的有效途徑，本文從彈性逆散射原理出發，首先推導了非均質HTI介質中縱波散射系數方程，并通過裂縫各向異性參數的引入，提出了一種新穎的AEI參數化方法，奠定了非均質HTI介質AEI反演的理論基礎。在此基礎上，發展了貝葉斯框架下柯西稀疏正則化及平滑模型正則化約束的HTI介質AEI反演方法，實現了裂縫弱度參數的穩定估算。模型和實際資料處理表明，該方法可從方位疊前地震資料中穩定估算裂縫弱度參數，為非均質HTI介質的各向異性預測提供了一種高可靠性的地震反演方法。

[1] Thomsen L.Understanding Seismic Anisotropy in Exploration and Exploitation.SEG 2010 Distinguished Instructor Short Course，2002.

[2] Tsvankin L and Grechka V.Seismology of Azimuthally Anisotropic Media and Seismic Fracture Characteri-zation.SEG Publication，USA，2011.

[3] Rüger A.P-wave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry.Geophysics，1997，62(3)：713-722.

[4] Rüger A.Variation of P-wave reflectivity with offset and azimuth in anisotropic media.Geophysics，1998，63(3)：935-947.

[5] Liu E and Martinez A.Seismic Fracture Characterization.EAGE Publication，Netherlands，2012.

[6] 陳懷震.基于巖石物理的裂縫型儲層疊前地震反演方法研究[學位論文].山東青島：中國石油大學(華東)，2015.

Chen Huaizhen.Study on Methodology of Pre-stack Seismic Inversion for Fractured Reservoirs Based on Rock Physics [D].China University of Petroleum (Huadong)，Qingdao，Shandong，2015.

[7] Hsu C J and Schoenberg M.Elastic waves through a simulated fractured medium.Geophysics，1993，58(7)：964-977.

[8] Bakulin A，Grechka V and Tsvankin I.Estimation of fracture parameters from reflection seismic data - Part Ⅰ：HTI model due to a single fracture set.Geophy-sics，2000，65(6)：1788-1802.

[9] Thomsen L.Weak elastic anisotropy.Geophysics，1986，51(10)：1954-1966.

[10] Tsvankin L.P-wave signatures and notation for transversely isotropic media：an overview.Geophysics，1996，61(2)： 467-483.

[13] 印興耀，宗兆云，吳國忱.非均質介質孔隙流體參數地震散射波反演.中國科學：地球科學，2013，43(2)：1934-1942.

Yin Xingyao，Zong Zhaoyun and Wu Guochen.Seismic wave scattering inversion for fluid factor of heterogeneous media.Science China：Earth Sciences，2013，43(12)：1934-1942.

[14] Wu R S and Aki K.Scattering characteristics of elastic waves by an elastic heterogeneity.Geophysics，1985，50(4)：582-595.

[15] Weglein A B.Nonlinear inverse scattering for multiple attenuation.Society of Photo-Optical Instrumentation Engineers Conference Series，1993，158-160.

[16] Stolt R H，Weglein A B and AbdElhadi Y E.Migration and inversion of seismic data.Geophysics，1985，50(12)：2458-2472.

[17] Mora P，Sarwar A K M and Smith D L.Nonlinear two-dimensional elastic inversion of multioffset seismic data.Geophysics，1987，52(9)：1211-1228.

[18] Shaw R K and Sen M K.Born integral，stationary phase and linearized reflection coefficients in weak anisotropic media.Geophysical Journal International，2004，158(1)：225-238.

[19] Shaw R K and Sen M K.Use of AVOA data to estimate fluid indicator in a vertically fractured medium.Geophysics，2006，71(3)：C15-C24.

[20] Connolly P.Elastic impedance.The Leading Edge，1999，18(4)：438-452.

[21] Whitcombe D N.Elastic impedance normalization.Geo-physics，2002，67(1)：60-62.

[22] 馬勁風.地震勘探中廣義彈性阻抗的正反演.地球物理學報，2003，46(1)：118-124.

Ma Jinfeng.Forward modeling and inversion method of generalized elastic impedance in seismic exploration.Chinese Journal of Geophysics，2003，46(1)：118-124.

[23] 李愛山，印興耀，陸娜等.兩個角度彈性阻抗反演在中深層含氣儲層預測中的應用.石油地球物理勘探，2009，44(1)：87-92.

Li Aishan，Yin Xingyao，Lu Na et al.Application of elastic impedance inversion with two angle stack ga-thers to predict gas-bearing reservoir of mid-deep layer.OGP，2009，44(1)：87-92.

[24] 張豐麒，金之鈞，盛秀杰等.一種穩健的彈性阻抗反演方法.石油地球物理勘探，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.

[25] Martins J L.Elastic impedance in weakly anisotropic media.Geophysics，2006，71(3)：2092-2096.

[26] 吳國忱，趙小龍，羅輯等.基于擾動彈性阻抗的裂縫參數反演方法.石油地球物理勘探，2017，52(2)：340-349.

Wu Guochen，Zhao Xiaolong，Luo Ji et al.Fracture parameter inversion based on perturbation elastic impe-dance.OGP，2017，52(2)：340-349.

[27] 羅輯，吳國忱，宗兆云等.基于方位彈性阻抗反演的裂縫型儲層流體檢測方法.石油地球物理勘探，2015，50(6)：1154-1165.

Luo Ji，Wu Guochen，Zong Zhaoyun et al.Fluid identification in fractured reservoirs based on azimuthal elastic impedance inversion.OGP，2015，50(6)：1154-1165.

[28] 陳懷震，印興耀，張金強等.基于方位各向異性彈性阻抗的裂縫巖石物理參數反演方法研究.地球物理學報，2014，57(10)：3431-3441.

Chen Huaizhen，Yin Xingyao，Zhang Jinqiang et al.Seismic inversion for fracture rock physics parameters using azimuthally anisotropic elastic impedance.Chinese Journal of Geophysics,2014，57(10)：3431-3441.

[29] Sacchi M D and Ulrych T J.High-resolution velocity gathers and offset space reconstruction.Geophysics，1995，60(4)：1169-1177.

[30] Alemie W and Sacchi M D.High-resolution three-term AVO inversion by means of a Trivariate Cauchy pro-bability distribution.Geophysics，2011，76(3)：R43-R55.

[31] 宗兆云，印興耀，吳國忱.基于疊前地震縱橫波模量直接反演的流體檢測方法.地球物理學報，2012，55(1)：284-292.

Zong Zhaoyun，Yin Xingyao and Wu Guochen.Fluid identification method based on compressional and shear modulus direct inversion.Chinese Journal of Geophysics，2012，55(1)：284-292.

[32] 宗兆云，印興耀，張峰等.楊氏模量和泊松比反射系數近似方程及疊前地震反演.地球物理學報，2012，55(11)：3786-3794.

Zong Zhaoyun，Yin Xingyao，Zhang Feng et al.Reflection coefficient equation and prestack seismic inversion with Young’s modulus and Poisson ratio.Chinese Journal of Geophysics，2012，55(11)：3786-3794.

[33] Scales J A and Smith M L.Introductory Geophysical Inverse Theory.Samizdat Press，1994.

[34] Schoenberg M.Elastic wave behavior across linear slip interfaces.Journal of the Acoustical Society of America，1980，68(5)：1516-1521.

[35] Schoenberg M.Reflection of elastic waves from periodically stratified media with interfacial slip.Geophysical Prospecting，1983，31(2)：265-292.

[37] Buland A and More H.Bayesian linearized AVO inversion.Geophysics，2003，68(1)：185-198.

[38] 張廣智，陳懷震，王琪等.基于碳酸鹽巖裂縫巖石物理模型的橫波速度和各向異性參數預測.地球物理學報，2013，56(5)：1707-1715.

Zhang Guanzhi，Chen Huaizhen，Wang Qi et al.Estimation of S-wave velocity and anisotropic parameters using fractured carbonate rock physics model.Chinese Journal of Geophysics，2013，56(5)：1707-1715.

附錄A非均質HTI介質縱波散射系數方程推導

式(10)中慢度矢量p和極化矢量t的下標m，n，i，j，k和l可表示為

(A-1)

式中δij和δkl均表示Kronecker符號。

根據Shaw等[19]給出的縱波慢度矢量p、極化矢量t及ηmn表達式，代入式(8)中，則非均質HTI介質的縱波散射系數方程可最終表示為

當介質為均質各向同性介質時，式(A-2)退化為各向同性介質的縱波散射系數方程[31]。

附錄B無量綱裂縫各向異性參數的引入

(B-1)

式中下標b表征法向裂縫弱度參數的均值屬性。對式(B-1)左右兩側進行積分，可得

(B-2)

(B-3)

式中m和n均為兩個未知量。將式(B-3)代入式(B-2)，可得

(B-4)

同樣，將式(B-4)代入式(B-3)，并取C=1，則

(B-5)

顯然，對裂縫弱度參數ΔT作類似變形處理，即得

(B-6)

*山東省青島市經濟技術開發區長江西路66號中國石油大學(華東)地球科學與技術學院，266580。 Email：panxinpeng1990@gmail.com

本文于2016年12月26日收到，最終修改稿于2017年10月6日收到。

本項研究受國家自然科學基金項目(41674130)、國家重點基礎研究發展計劃項目(2013CB228604，2014CB239201)、國家油氣重大專項(2016ZX05027004-001、2016ZX05002005-09HZ)、中央高?；究蒲袠I務費專項資金課題(15CX08002A、17CX06040)等聯合資助。

1000-7210(2017)06-1226-10

潘新朋，張廣智，印興耀.非均質HTI介質裂縫弱度參數地震散射反演.石油地球物理勘探，2017,52(6):1226-1235.

P631

A

10.13810/j.cnki.issn.1000-7210.2017.06.013

(本文編輯：朱漢東)

潘新朋 博士研究生，1990年生； 2013年獲中國石油大學(華東)地球物理學專業理學學士學位； 2016年獲中國石油大學(華東)地球物理學專業理學碩士學位； 現在中國石油大學(華東)地球科學與技術學院攻讀地質資源與地質工程專業博士學位，主要從事裂縫儲層巖石物理、儲層預測與流體識別等方面的學習和研究。