基于Ruger 近似式預測裂縫型儲層的可行性研究

吳姍姍，高 剛，桂志先，王 鵬，李成毅

(1.長江大學油氣資源與勘探技術教育部重點實驗室，湖北 武漢 430100；2.長江大學地球物理與石油資源學院，湖北 武漢 430100；3.中國石油集團川慶鉆探工程有限公司地球物理勘探公司，四川 成都 610213)

裂縫型儲層的預測方法根據資料來源通常分為兩類：一類是通過測井資料與巖心露頭資料定量識別儲層的裂縫，該方法預測結果精度較高，但其結果難以外推至井間[1]；另一類是通過地震資料定性預測儲層的裂縫，由于地震資料有橫向分辨率高、勘探成本較低、數據分布范圍廣等特點，因而運用地震資料檢測裂縫的方法廣泛應用于裂縫型儲層的預測[2]。目前運用地震資料檢測裂縫的方法主要有：①疊后地震屬性分析法，即利用疊后資料分析其地震屬性，如相干屬性[3]、曲率屬性[4]、螞蟻體屬性[5]等；② 轉換波裂縫檢測法，利用橫波分裂產生的快慢橫波時差來反映裂縫發育的主方向和裂縫密度，常用的方法有相對時差梯度法[6]和層剝離法[7]；③疊前縱波方位各向異性檢測法，利用縱波各向異性進行裂縫檢測的方法有動校正(Normal Moveout Correction，簡稱NMO)速度方位變化裂縫檢測法、正交地震測線縱波時差裂縫檢測法、縱波方位AVO(Amplitude Variation with Offset and Azimuth，簡稱AVOZ/AVAZ)和縱波阻抗隨方位角變化(Impedance Variation with Azimuth，簡稱IPVA)裂縫檢測法[8]。

利用疊前縱波資料的方位各向異性特征來識別裂縫型儲層是目前應用最為廣泛的方法之一。該技術的研究始于20 世紀70 年代，由于受到當時勘探條件的限制，直到20 世紀90 年代，I.Tsvankin[9]和A.Ruger[10]分別推導了HTI(Horizontal Transverse Isotropy，簡稱HTI)介質中純模式波的動校正速度隨方位角變化(Velocity Variation with Azimuth，簡稱VVA)和縱波的反射系數隨入射角和方位角變化(AVOZ/AVAZ)的近似公式之后，該技術才引起了廣泛的關注。隨后S.Malick 等[11]、范國章等[12]、朱成宏等[13]、朱兆林等[14]對地震波方位AVO 特性進行了全面的分析，建立了地震波方位AVO 的變化特征與裂縫走向、傾向以及發育強度的關系，該方法在實際裂縫檢測應用中也取得了很好的效果[15-19]。該方法在具體應用時僅考慮了各向異性參數對地震響應特征的影響，并未考慮巖性參數對地震響應特征的影響。然而，在實際資料的應用過程中發現，在儲層與蓋層的彈性參數存在較為嚴重疊置現象的情況下，常規的彈性參數預測儲層難度較大，雖然裂縫參數預測縱波方位各向異性的相關技術方法可以提高儲層的預測精度，但是存在一定的多解性。因此，筆者以Ruger 近似式為基礎建立模型，將模型中各向異性參數和巖性參數設置為概率密度分布函數，通過Monte Carlo 隨機方法進行疊前AVAZ正演模擬擬合橢圓，進而分析各參數的隨機分布對擬合的橢圓信息即裂縫信息的影響。

1 疊前AVAZ 技術檢測裂縫原理

1.1 常用裂縫模型

在眾多的裂縫巖石物理模型中，目前最為常用的模型是Hudson 的薄硬幣模型和Schoenberg 的線性滑動模型。

Hudson 模型主要是描述彈性固體中包含薄硬幣形狀的橢球裂縫的等效介質模型，該模型主要通過裂縫密度e和裂縫的縱橫比χ來表征裂縫系統，裂縫的密度是采用裂縫體積密度的定義方式：

式中a是裂縫的半徑；N/V是單位體積內裂縫的數量；φ是裂縫誘導的孔隙率；χ是裂縫的縱橫比。

Schoenberg 線性滑動模型是基于Backus 細層層狀介質模型提出的。該模型在巖石物理模型的研究中，適用于具有線性連續邊界的、充滿弱強度(充填物彈性模量小)介質的平行層模型。C.J.Hsu 等[20]提出了兩個無量綱的參數，用來描述由裂縫引起的巖石剛度變化特征。

式中ΔN和ΔT分別代表由裂縫引起的垂直于裂縫面方向和平行于裂縫面方向的彈性參數差值，它們的變化范圍為0～1；EN和ET為非負的、無量綱的與裂縫有關的柔度參數；ZN為垂向柔度；ZT為平行于裂縫面的柔度；μ、λ為拉梅參數。

兩種模型參數之間具有一定的關聯性，ΔN和ΔT與裂縫參數之間的關系見式(4)和式(5)。

式中g=μ/(λ+2μ)；K′和μ′是縫隙中充填物的體積模量和剪切模量。

A.Bakulin 等[21]給出了HTI 介質中各向異性參數(ε(V)、γ、δ(V))與Schoenberg 裂縫模型巖石物理參數(ΔN和ΔT)之間的關系，如式(6)—式(8)所示：

將式(4)—式(5)代入式(6)—式(8)可以得到HTI介質中各向異性參數(ε(V)、γ、δ(V))與Hudson 模型所定義的裂縫密度之間的關系，如式(9)—式(11)所示，揭示了HTI 介質中各向異性參數(ε(V)、γ、δ(V))與裂縫的密度之間的關系，可以分析L.Thomsen[22]各向異性參數與裂縫密度及裂縫中流體之間的關系。

1.2 HTI 介質中疊前AVAZ 技術檢測裂縫原理

A.Ruger 等[10]對HTI 介質中地震波傳播理論(圖1)進行了深入研究后，在弱各向異性的假設條件下，

給出了HTI 介質中縱波的反射系數隨入射角和方位角變化的近似公式，即：

式中Z為縱波波阻抗；β為縱波波速；α為橫波波速；(ε(V)、δ(V)為不同于L.Thomsen 裂隙參數，后續為了方便，參數直接寫為ε、δ)，ε、γ分別表示縱波和橫波各向異性；δ表示縱向和橫向各向異性；G=ρβ2表示切向模量；i和φ分別為入射方向和測線方向與裂縫對稱軸的夾角；符號“-”和“Δ”分別表示上下界面參數的算術平均值、差值。

HTI 介質的各向異性參數的表達式：

式中cij為HTI 介質彈性常數。

圖1 HTI 介質中縱波傳播示意圖Fig.1 Sketch map of P wave propagation in HIT medium

入射角較小時，各向異性介質中縱波反射振幅與裂縫參數之間的關系可以近似表示為：

式中P為縱波垂直入射時的反射振幅；Giso為反射振幅隨偏移距的變化率，即各向同性梯度；Gani為振幅隨偏移距和方位角的變化率，即各向異性梯度，指示了介質受各向異性影響的部分。

在固定入射角的情況下，式(16)可以進一步簡化為：

式中A為與炮檢距有關的偏置因子；B為與炮檢距和裂縫特征有關的調制因子；θ是炮檢距方向與裂縫走向的夾角。

式(17)反射系數R隨著角度θ的變化可用橢圓狀的圖形來表示(圖2)。

參數A與炮檢距有關，反映了巖性變化所引起的振幅的變化量；參數B相當于振幅隨炮檢距方位改變的變化量，反映了地層裂縫對反射振幅強度的影響，其大小決定了儲層裂縫的發育程度，可視為各向異性因子。因此，可以把B/A即橢圓的扁率作為裂縫密度的相對量度，B/A的值越大，說明裂縫發育越好，反之，裂縫發育越差，而對于各向同性介質，A值就是均勻介質下的反射振幅，B值為零。

圖2 地震反射振幅隨方位角變化Fig.2 Variation of seismic reflection amplitude with azimuth

2 Monte Carlo 隨機正演模擬

疊前AVAZ 技術檢測裂縫法主要是通過地震資料進行橢圓擬合，然后根據橢圓的信息對裂縫進行預測和評價。在實際應用中，當橫向巖性變化時，通常是利用平均振幅A的信息，沒有考慮橫向巖性變化對擬合橢圓信息的影響，這在一定程度上造成了多解性。為了了解橫向巖性參數的變化對橢圓擬合信息的影響，對儲層物性較好，以孔-縫型、裂縫型儲層為主的研究區域的測井資料進行統計，灰巖儲層的縱波速度為3 000～5 200 m/s，密度為2.40～2.75 g/cm3，蓋層泥巖縱波速度為2 600～4 500 m/s，密度為2.1～2.5 g/cm3，據此建立含EDA 介質的雙層模型(圖3)。模型參數設置見表1，入射角為30°，采用Monte Carlo 隨機方法產生滿足條件的隨機模型，利用最小二乘法進行橢圓擬合，然后分析Ruger公式(式(12))中各個參數變化對擬合橢圓信息的影響，進而說明該參數對裂縫評價的影響。

圖3 含EDA 介質的雙層介質模型Fig.3 Two-layer model with EDA medium

2.1 各向異性參數對疊前AVAZ 地震響應影響

首先，通過模型分析各向異性參數單獨變化對疊前AVAZ 地震響應特征的影響。L.Thomsen[23]通過實驗測定了一系列不同巖性的樣品，根據其測定的不同巖性的各向異性參數值進行統計，發現各向異性參數值基本都小于0.2，所以本文模型各向異性參數分別設置為均值0.05、0.10、0.20 的正態分布，其中標準差設置為均值的10%，為了保證擬合出的橢圓扁率變化是由各向異性參數變化所引起，將模型速度和密度設置為固定的平均值(表1)，利用Ruger 近似式進行隨機正演模擬，分別作出不同各向異性參數值所對應的橢圓扁率B/A與各向異性因子B散點圖(圖4)。

表1 雙層介質模型參數Table 1 Parameters of a two-layer medium model

圖4 不同各向異性參數的B/A(橢圓扁率)與B(各向異性因子)散點圖Fig.4 The scatter plots of B/A and B with different anisotropic parameters

圖4a—圖4c 分別對應各向異性參數ε(圖4a)、γ(圖4b)、δ(圖4c)為不同均值的正態分布時的橢圓的扁率B/A與各向異性因子B散點圖，其中紅色、紫色及藍色分別為模型的各向異性參數值設置為均值0.05、0.10 和0.20 的正態分布結果。由圖4 可以看出橢圓的扁率B/A與各向異性因子B的值隨各向異性參數值增大而增大，說明各向異性越強，其各向異性因子B與橢圓的扁率B/A值越大，證明了利用縱波方位的各向異性特征檢測裂縫是可行的。對比圖4a—圖4c 可以看到，各向異性參數在同一均值的條件下，參數γ(圖4b)對應的橢圓扁率B/A與各向異性因子B的值最大，其次是參數δ(圖4c)，參數ε(圖4a)最小，因此，各向異性參數γ對橢圓扁率B/A與各向異性因子B影響最大、δ次之、ε最小。

2.2 非各向異性參數對疊前AVAZ 地震響應影響

為了進一步分析非各向異性參數對疊前AVAZ地震響應特征的影響，以各向異性參數γ正態分布為例(各向異性參數γ對疊前AVAZ 響應特征影響最大，模型參數γ設置為均值0.1、標準差0.01 的正態分布)，然后分別將模型中所研究的巖性參數根據表2 設置為正態分布，標準差按照模型4 取值，最后通過Monte Carlo 隨機正演模擬得到對應橢圓扁率B/A與各向異性因子B散點圖(圖5)。

圖5a 和圖5d 分別是上層和下層的密度設置為正態分布(表2 中模型4)，縱波和橫波速度設置為均值時所對應的橢圓扁率B/A與各向異性因子B散點圖，可以看出在密度標準差接近均值5%的情況下，其橢圓扁率B/A與各向異性因子B散點圖形態與圖4b 中γ為均值0.1 模型較為接近，說明密度參數不是影響橢圓信息的主控因素。圖5b 和圖5e 分別為上層介質和下層介質的橫波速度設置為正態分布(表2 中模型4)，縱波速度和密度設置為均值時所對應的橢圓扁率B/A與各向異性因子B散點圖，當橫波速度的標準差接近均值的10%時，其圖形的散點離散程度較密度正態分布時大，與圖4b 中γ為均值0.1 的模型結果差異較大，可以看出橫波速度隨機變化對AVAZ 響應特征的影響較大。圖5c 和圖5f 分別為上層介質和下層介質的縱波速度設置為正態分布(表2 中模型4)，橫波速度和密度設置為均值時所對應的橢圓扁率B/A與各向異性因子B散點圖，綜合來看，當縱波速度所取的標準差接近均值的10%時，其圖形的散點擴散程度最大，與圖4b 中γ為均值0.1 的模型結果差異很大，說明縱波速度的隨機變化對AVAZ 響應特征的影響最大。

表2 不同標準差的巖性參數表Table 2 Lithology parameters with different standard deviations

圖5 不同非各向異性參數B/A(橢圓扁率)與B(各向異性因子)散點圖(γ 隨機變化、ε=0、δ=0)Fig.5 The scatter plots of B/A and B with different non-anisotropic parameters(γ random variation,ε=0,δ=0))

對式(12)進行理論分析，可以看出公式中的密度參數對橢圓的擬合影響較小，速度參數影響較大。式(12)的理論分析結果與上述模型分析結果相符合，因此，模型速度與密度的不確定性會對模型模擬結果產生影響，其中縱波速度的影響最大，橫波速度的影響次之，密度的影響較小。

由于在實際地層中速度和密度并不是單獨變化，為了進一步探究非各向異性參數的隨機變化對隨機正演的影響程度，在各向異性參數ε(圖6a、圖7a、圖8a)、γ(圖6b、圖7b、圖8b)、δ(圖6c、圖7c、圖8c)分別設置為均值為0.1，標準差為0.01 的正態分布的基礎上，將地層的速度和密度按照表2(模型1—模型4)設置為不同標準差的正態分布，采用Monte Carlo 隨機方法產生滿足條件的隨機模型，利用Ruger 近似式正演其疊前AVAZ 地震響應，生成橢圓扁率B/A與各向異性因子B散點圖(圖6—圖8)。其中，圖6 為上層介質的速度與密度按照模型1—模型4 設置為不同標準差的正態分布，下層介質的速度與密度設置為恒定的均值；圖7 為下層介質的速度與密度按照模型1—模型4 設置為不同標準差的正態分布，上層介質的速度與密度設置為恒定的均值；圖8 為上、下兩層介質的速度與密度均按照模型1—模型4 設置為不同標準差的正態分布。

圖6 不同標準差的速度與密度模型B/A(橢圓扁率)與B(各向異性因子)散點圖(上層介質速度與密度均隨機變化)Fig.6 B/A and B scatter plots for velocity and density models with different standard deviations (upper layer velocity and density are set to a normal distribution)

圖7 不同標準差的速度與密度模型B/A(橢圓扁率)與B(各向異性因子)散點圖(下層介質速度與密度均隨機變化)Fig.7 B/A and B scatter plots for velocity and density models with different standard deviations(The velocity and density settings of the upper layer are randomly distributed)

圖8 不同標準差的速度與密度模型B/A(橢圓扁率)與B(各向異性因子)散點圖(兩層介質速度與密度均隨機變化)Fig.8 B/A and B scatter plots for velocity and density models with different standard deviations (both the velocity and density of two-layer medium change randomly)

在圖6—圖8 中均可看出：速度與密度設置為均值時(模型1)，其橢圓扁率的變化完全是由各向異性參數的隨機變化所引起，由于模型2—模型4 的速度與密度設置為正態分布，其擬合的橢圓信息的變化是由巖性參數和各向異性參數共同變化所引起的結果。分別將模型2—模型4 與模型1 對比分析：速度與密度在較小的范圍內(模型2)隨機取值時，模擬出的橢圓信息可以很好地反映裂縫的信息；速度與密度的標準差分別為均值的5%和2.5%時(模型3)，模擬出的散點圖基本可以反映裂縫的信息；但當速度與密度的標準差分別為均值的 10%和 5%時 (模型4)，橢圓信息很大一部分受到速度和密度隨機波動的影響，其散點分布較為離散。

對比圖6–圖8 可以看出，兩層介質(圖8)的速度與密度均為正態分布，比僅上層介質(圖6)或僅下層介質(圖7)的速度與密度為正態分布時所模擬出的橢圓扁率B/A與各向異性因子B散點圖更離散，由于樣本值滿足正態分布，其值大部分仍然分布在均值附近，所以其散點密集區基本一致。圖6—圖8 中間一列的b 圖為各向異性參數γ設置為概率密度分布，其B/A與B散點分布范圍最廣；圖6—圖8 左側一列的a 圖為各向異性參數ε設置為概率密度分布，其橢圓扁率B/A與各向異性因子B散點圖分布范圍最小；圖6—圖8 右側一列c 圖為各向異性參數δ設置為概率密度分布時，其橢圓扁率B/A與各向異性因子B散點圖分布范圍介于a 圖和b 圖之間。分析得出如下結論：

a.兩層介質的速度與密度均設置為正態分布，較單層介質(上層介質或下層介質)的速度與密度正態分布對橢圓扁率B/A與各向異性因子B散點圖影響更大。

b.無論是儲層還是蓋層，由橫向巖性變化大而引起速度與密度的隨機變化時，橢圓的扁率B/A與各向異性因子B的變化范圍均會變大。當橫向巖性變化過大即速度與密度隨機波動過大(例如模型 4 的速度隨機變化范圍是其平均值的10%，密度隨機變化范圍是其平均值的 5%)，會無法準確識別出由各向異性所引起的變化。在不同巖性的裂縫性儲層中，地層的橫向巖性變化較小時即速度與密度隨機變化范圍在較小的范圍內(標準差小于均值的10%)才能較好地利用該技術評價裂縫。

3 結論

a.以含EDA 介質的雙層模型為基礎，將Ruger近似式中的各向異性參數和巖性參數設置為正態分布，通過Monte Carlo 隨機方法進行疊前AVAZ 正演模擬擬合橢圓信息，分別分析了各向異性參數和巖性參數的隨機變化對擬合的橢圓的扁率B/A和各向異性因子B產生影響。模擬實驗結果得知，在利用AVAZ 技術評價HTI 介質中裂縫發育程度時，不僅各向異性參數ε、γ和δ會對AVAZ 響應擬合的橢圓扁率B/A和各向異性因子B產生影響，而且非各向異性參數如地層的速度和密度(尤其是地層的縱波速度)的隨機變化也會對擬合的橢圓扁率B/A和各向異性因子B造成一定的影響。

b.各向異性參數γ對橢圓扁率B/A和各向異性因子B的影響最大，δ次之，ε最小。非各向異性參數對橢圓扁率B/A和各向異性因子B的影響程度依次為，縱波速度最大、橫波速度次之，密度最小。因此，縱波速度的不確定性是造成AVAZ 的響應特征分布的主控因素，橫波速度的影響也較大，密度的影響較小。當速度和密度的隨機波動范圍過大時，與理想的僅考慮各向異性參數變化的模型結果差別較大，所以當速度與密度在較小的范圍內隨機取值時，橢圓扁率B/A和各向異性因子B才可以較好地反映地層各向異性的強弱即裂縫的密度大小。

c.應用疊前縱波方位各向異性檢測儲層裂縫時，應當充分考慮地層速度和密度橫向的變化對AVAZ 響應特征分布的影響。在具體的實際資料分析中，為了減少地層速度、密度參數對預測結果的影響，在疊前三參數反演得到空間的縱橫波速度、密度的基礎上，下一步可研究如何利用空間的縱橫波速度、密度來提高橢圓扁率的計算精度，進而提高該參數預測裂縫信息的準確性。