流體飽和正交裂縫孔隙型介質數值模擬及波場特征分析

王鈺鵬, 李小斌, 閻建國

(成都理工大學 地球物理學院,成都 610059)

0 引言

裂縫是油氣重要的運移通道和儲存空間，因此裂縫發(fā)育的地層是良好的油氣儲集層。對裂縫介質各向異性的研究，主要是通過等效介質理論建立裂縫參數與地震各向異性之間的聯系。相對發(fā)展成熟的傳統等效介質理論有Hudson擾動理論[1]、Thomsen 等徑孔隙理論[2]、線性滑動理論[3]和Chapman多尺度裂縫等效介質[4]。

近年來，許多學者充分考慮介質的裂縫、孔隙以及之間的流體交換作用，提出了一系列的等效介質模型。Parra[5]結合BISQ噴流機制和Thomsen模型；分析在不同方位角和入射角傳播的地震波的速度頻散和衰減特性;Gurevich[6]將線性滑動理論和Gassmann各向異性理論相結合，推導了垂直裂隙和孔隙液壓連通的等效模型；杜啟振等[7]利用線性滑動理論、Biot孔隙介質理論和連續(xù)介質滲流模型建立了裂縫誘導各向異性雙孔隙等效模型，推導了速度應力方程，并使用有限差分方法進行了數值模擬；Bakulin等[8]依據線性滑動理論，給出了正交介質的彈性系數計算方法；張顯文等[9]利用BISQ機制推導了雙相正交各向異性介質的頻散方程，研究了地震波的衰減頻散方位特性；高煒等[10]將Biot雙相介質理論和Gurevich模型相結合，推導了裂縫誘導雙相介質模型的彈性波傳播方程；劉財等[11]使用偽譜法對高煒建立的模型進行數值模擬，研究其波場特征；張繁昌等[12]將裂縫看作高孔隙薄層周期層狀模型，結合Biot方程，推導了裂隙介質的剛度系數與頻率的關系，分析孔隙和裂縫之間的不同類型的擠噴流影響。

為了研究正交裂縫，孔隙并存時，裂縫參數對波場傳播特征的影響，筆者參照了Gurevich的思路，首先利用線性滑動理論，建立兩組相互正交的裂縫，然后通過各向異性Gassmann方程，引入孔隙背景，最后使用Biot雙相介質理論彈性波傳播方程，用有限差分方法進行數值模擬。通過分析波場特征，進一步認識裂縫性儲層的地震波傳播規(guī)律。筆者所建立的這種模型旨在模擬和分析“裂縫、孔隙、流體共存”介質中的波場特征，為地震儲層預測及方法研究提供一些參考。

1 基本原理

1.1 孔隙介質

本次研究所建立的模型如圖1所示，含兩組相互正交的裂隙，一組裂隙法向平行與x軸的豎直裂隙，另一組法向平行與z軸的水平裂隙。在實際儲層中，豎直裂隙常常由構造運動造成，水平裂隙提供了流體流動和儲集空間。模型的背景巖石為各向同性高孔隙性介質。

圖1 正交裂隙孔隙型介質Fig.1 Orthogonal fractured porous medium

在各向同性介質中，引入均勻分布的孔隙，結果為含孔隙的各向同性介質。使用以下經驗公式來計算引入孔隙后的彈性模量[13]，

K=(1-α0)Kg

(1)

μ=(1-α0)μg

(2)

(3)

式中：k、μ分別為含孔隙巖石體積模量和剪切模量;kg、μg分別為背景巖石顆粒的體積模量和剪切模量;φp孔隙的孔隙度。

1.2 正交裂縫

按照Schoenberg理論，整個介質的彈性柔度張量表示為式(4)。

(4)

(5)

另一組旋轉不變性裂隙法向平行于z軸時，設n3=(0，0，1)，則柔度張量可以寫成式(6)形式。

(6)

因此，當介質中同時含有一組法向平行于x軸和一組法向平行于z軸的裂縫時，介質的柔度張量可以表示為式(7)。

S0=Sb+Sf1+Sf2

(7)

對式(7)求逆可得正交裂縫介質的彈性系數。

C11=(λ+2μ+4μ(λ+μ)ZN2)/[(λ+2μ)·

(ZN1+ZN2)+4μ(λ+μ)ZN1ZN2+1]

C12=(λ+2λμZN2)/[(λ+2μ)(ZN1+ZN2)+

4μ(λ+μ)ZN1ZN2+1]

C13=λ/[(λ+2μ)(ZN1+ZN2)+

4μ(λ+μ)ZN1ZN2+1]

C22=[λ+2μ+4μ(λ+μ)(ZN1+ZN2)+

(8μ3+12λμ2)ZN1ZN2]/[(λ+2μ)·

(ZN1+ZN2)+4μ(λ+μ)ZN1ZN2+1]

C23=(λ+2λμZN1)/[(λ+2μ)(ZN1+ZN2)+

4μ(λ+μ)ZN1ZN2+1]

C33=(λ+2μ+4μ(λ+μ)ZN1)/[(λ+

2μ)(ZN1+ZN2)+4μ(λ+μ)ZN1ZN2+1]

(8)

式中:λ、μ為含孔隙各向同性背景介質的拉梅系數和剪切模量。

1.3 孔隙流體

研究飽和流體裂隙孔隙介質中的波場傳播，可以利用各向異性Gassmann方程[14]引入孔隙流體。當介質由單一各向同性顆粒組成的各項異性骨架時，流體飽和彈性模量和干巖石骨架彈性系數之間的關系可以用Gurevich改的各向異性Gassmann方程來表示[15]:

(9)

其中，

式中:Kg為組成巖石的各向同性顆粒的體積模量;Kf為孔隙流體的體積模量;φ是總孔隙度。

1.4 Biot模型

1.4.1 本構方程

Biot理論基于一定的假設描述了流體飽和孔隙介質中的彈性波傳播。根據Biot的描述，本構方程為式(10)。

(10)

其中：ekl為固相應變；ε為液相應變；Qj表示固液之間的體積變化耦合關系；R為孔隙流體的彈性參數；σij、s分別為固、液相的應力，并且孔隙度φ及流體壓力p滿足,

s=-φp

(11)

1.4.2 運動微分方程

根據Biot理論，不考慮外力作用時，飽和雙相介質的運動微分方程為式(12)。

(12)

式中:ρ11+ρ12=(1-φ)ρs,ρ12+ρ22=φρf，ρ12=-ρa

其中：B為耗散系數;η為流體的粘滯系數;kij為各個方向的滲透率;ρs為固相骨架密度;ρf為液相流體密度;ρa為兩者之間的耦合密度。

1.4.3 幾何方程

若固相位移為u，液相位移為U，在雙相介質中，應力與位移的關系有

(13)

1.4.4 雙相介質的一階速度應力方程

將幾何方程(13)帶入到本構方程(10)，再對時間求一階倒，得到速度表示應力的方程為式(14)。

(14)

將運動微分方程(12)化簡，得到應力表示速度的方程為式(15)。

(15)

2 數值模擬

2.1 模型設置和網格參數

我們主要研究波場隨著裂縫參數的變化，因此除裂縫外的參數是不變的，參數值見表1。

表1 介質參數表(不包括裂縫參數)

使用有限差分方法正演 ，震源使用雷克子波，函數如下，

e-λ[(x-x0)2+(z-z0)2]

(16)

密度單位kg·m-3，耦合參數Q1、Q3、流體參數R單位為109kg·m-1，Kf單位為Gpa，速度Vp和Vs單位為m/s

式中:f0為震源中心頻率;t0=1/f0為子波延遲時間;(x0,z0)為震源中心位置;λ為震源作用的集中系數。

設計一個單層均勻的正交裂縫雙相模型，模型網格數為300×300，網格大小為1 m×1 m。子波頻率為20 Hz，時間采樣間隔時0.1 ms。

2.2 不同裂隙參數的數值模擬

筆者使用裂隙參數為裂隙法向弱度EN和切向弱度ET。弱度參數的大小與裂隙尺度，填充物等參數相關。弱度與屈服度之間的關系為，

EN=ZN(λ+2μ),ET=ZTμ

(17)

式中:λ、μ為背景巖石介質的拉梅常數和剪切模量。

需要注意的是，文中的裂隙參數均為干裂隙參數，流體的表示是使用Gassmann方程來填充流體。在一些文章中提出EN為“0”表示流體飽和情況，實際上在高頻情況下(實驗室超聲波測量)才滿足這一條件，而本文通過Gassmann方程引入流體，表述的是低頻段(地震頻段)時候的飽和情況。若裂隙中沒有其他固體包含物，則裂隙的法向弱度和切向弱度之比僅僅與背景介質的縱橫波速度之比相關。

2.2.1 不同分量波場

按照表2的參數計算出四個分量的波場快照，如圖2所示。從圖2可以看出，雙相各向異性介質中有三種波，傳播速度最快的快qP波，呈現橢圓形，中間一圈形狀復雜，有三叉區(qū)現象的波是qSV波，最內一圈是傳播速度最慢的慢qP波。從流相與固相的對比中可以看出，流相的慢qP波比固相的慢qP波能量強，固相的qSV波比流相的qSV波能量強。

圖2 不同分量波場快照(0.04 s)Fig.2 Snapshots of wavefields of different components(0.04 s)(a)固相速度x分量;(b)固相z分量；(c)流相x分量;(d)流相z分量

表2 模型裂隙參數表

筆者主要是分析波場特征隨著裂縫參數的變化，考慮到各個分量的波形差異不大，后面只展示波場的固相速度z分量。

2.2.2 單組裂隙變化對波場特征的影響

在單組裂縫情況下，研究ET參數變化對波場特征的影響，裂縫參數的選擇如表3所示。計算出的波場快照如圖3所示。

表3 單組裂隙變化裂隙參數表

從圖3可以看出，隨著切向弱度ET的增加，快qP波在裂縫的兩側和裂縫走向上，即水平和豎直方向上波場位置沒有變化。在四個角處向內微微收斂，從跑道形態(tài)的橢圓轉向偏菱形的橢圓，變化不太明顯?？捎^察到qSV波的變化劇烈，從沒有三叉區(qū)現象到三叉區(qū)現象明顯。在介質的各個方向上波形縮小，速度降低。慢qP波形態(tài)變化不明顯，能量加強。

圖3 豎直裂隙切向弱度變化的波場快照(0.04 s)Fig.3 Snapshot of wavefields of vertical fracture tangential weakness change (0.04 s)(a)ET1=0.2分量;(b)ET1=0.4;(c)ET1=0.7;(d)ET1=1.0

2.2.3 單組裂隙變化對波場特征的影響

研究豎直裂縫模型波場隨著裂縫法向弱度參數的變化，設計裂縫參數如表4所示，計算出的波場快照如圖4所示。

圖4 豎直裂隙法向弱度變化的波場快照(0.04 s)Fig.4 Snapshot of wavefield of weakness change in normal direction of vertical fracture (0.04 s)(a)EN2=0.2;(b)EN2=0.4;(c)EN2=0.7;(d)EN2=1.0

表4 單組裂隙變化裂隙參數表

從圖4中可以看出，隨著的增加，qSV波變化不明顯，形狀呈現橢圓狀，存在不明顯的三叉區(qū)現象，快qP波在裂縫的兩側方向向內收斂，橢圓形態(tài)加重，長軸變長，短軸變短，扁平率增加，在裂縫的走向上沒有變化，傳播速度不變。慢qP波變化不明顯，在裂縫兩側有微弱的向內收斂。

2.2.4 正交裂隙變化對波場特征的影響

在正交裂縫情況下，研究波場快照隨著裂隙切向弱度的變化。兩組裂隙的參數值相同，具體參數見表5，計算出的波場快照如圖5所示。

表5 正交裂隙變化裂隙參數表

圖5 正交裂隙模型切向弱度變化的波場快照(0.04 s)Fig.5 Wavefield snapshot of tangent weakness change of orthogonal fracture model (0.04 s)(a)ET=0.1;(b)ET=0.3;(c)ET=0.6;(d)ET=1.0

從圖5中可以看出整體上，在ET較小的時候，三叉區(qū)現象明顯，位于四個坐標軸方向，如圖5(a)所示；隨著ET增加，三叉區(qū)現象減弱，然后再增強，此時的三叉區(qū)位于坐標軸的四個近似45°角處?？靟P波在水平和豎直方向上波場位置沒有變化，在四個角處向內微微收斂。qSV波變化劇烈，各方向的速度降低，波形縮小；三叉區(qū)現象的位置改變，從四個正方向到四個斜方向。慢qP波形態(tài)變化不明顯，能量加強。

2.2.5 正交裂隙變化對波場特征的影響

在正交裂縫情況下，研究波場快照隨著裂隙法向弱度的變化。裂隙參數見表6，計算出的波場快照如圖6所示。

圖6 正交裂隙模型法向弱度變化的波場快照(0.04 s)Fig.6 Snapshot of wavefield of normal weakness change of orthogonal fracture model (0.04 s)(a)EN=0.2;(b)EN=0.5;(c)EN=0.8;(d)EN=1.0

表6 正交裂隙變化裂隙參數表

從圖6中觀察到，快qP波在形態(tài)上都近乎橢圓，隨增大，橢圓縮小，形成較為標準的圓形。快qP波各個方向的速度降低。qSV波在水平和豎直方向上波場位置沒有變化。在四個角處的三叉區(qū)現象消失，波形向內靠。形態(tài)上從不規(guī)則到橢圓狀態(tài)。慢qP波形態(tài)變化不明顯，振幅增強。

從圖3到圖6，慢qP波的形態(tài)變化都不大，說明裂隙參數對慢qP波的影響較小。

2.2.6 裂隙密度對波場特征的影響

在上面的一些規(guī)律分析中，裂縫參數比較抽象，且各別極端情況并不會出現。筆者只討論一種具體的情況，波場隨著裂隙密度的變化的規(guī)律。從Hudson理論中，當裂縫為橢球體裂縫時，不含包含物的干裂縫弱度計算公式為：

(18)

按照本文設計的參數，計算出裂隙弱度參數隨著裂隙密度的變化，如圖7所示。法向弱度和切向弱度隨著裂縫密度線性增加。法向弱度相對于切向弱度隨裂隙密度的變化更敏感。這表明在做裂縫預測工作時，可以更多的關注法向弱度這個參數。

圖7 弱度隨裂縫密度變化曲線Fig.7 The curve of weakness-fracture density

按照不同的裂隙密度，計算出波場快照，單組裂隙模型如圖8，正交裂隙如圖9所示。

圖8 單組裂隙波場隨密度變化快照(0.04 s)Fig.8 Snapshot of the wavefield of a single group of fractures changing with density (0.04 s)(a)e=0.01裂縫密度;(b)e=0.05;(c)e=0.10;(d)e=0.15

圖9 正交裂隙波場隨密度變化快照(0.04 s)Fig.9 Snapshot of orthogonal fracture wavefield changing with density (0.04s)(a)e=0.01裂縫密度;(b)e=0.05;(c)e=0.10;(d)e=0.15

從模擬結果來看，單組裂縫時，快qP波在裂縫兩邊向內收斂，從圓形變?yōu)闄E圓；qSV波，各方向的速度降低，波形縮小，裂隙密度增大到一定程度，出現較弱的三叉區(qū)現象；正交裂縫時，qSV波和快qP波的形態(tài)近乎圓形，隨著裂隙密度增加，圓形縮小。

3 結論與展望

筆者使用線性滑動理論和各向異性Gassmann方程，構建飽和的正交裂縫孔隙型介質，然后使用高階有限差分對我們構建的正交模型進行數值模擬，并取得了以下幾點認識：

1)裂縫參數的變化對慢qP波的形態(tài)影響不大，但裂縫參數值增大，慢qP波能量增強。

2)EN變化對快qP波影響大。對qSV波的影響小；ET變化對qSV波影響大，對于快qP波而言，主要是影響除切向法法向以外的速度。

3)在單組裂隙中，隨著ET越大，三叉區(qū)現象約明顯。隨著EN的增加，快qP波的橢圓形態(tài)加重。

4)在兩組裂隙參數相同的正交裂隙介質中，隨著ET增加，qSV三叉區(qū)現象從四個坐標軸方向逐漸消失，然后再出現在坐標軸的四個近似45°角處。EN增加，快qP波縮小，形態(tài)趨近于標準的圓形。

在實際的地震勘探中，主要還是應用縱波勘探。通過本文的研究，縱波主要受到裂隙法向弱度的影響。若不考慮裂隙中含有其余固相顆粒，則對裂隙密度更為敏感，也就是對裂隙發(fā)育更敏感。這表明對于地下裂隙的預測，我們可以更多地關注如何從縱波地震數據中反演出裂隙法向弱度參數。