基于三維孔隙網絡模型的縱波頻散衰減特征分析

魏樂樂， 甘利燈*， 熊繁升， 孫衛濤， 丁騫， 楊昊

1 中國石油勘探開發研究院， 北京 100083 2 清華大學周培源應用數學研究中心， 北京 100084

0 引言

地震波在通過孔隙介質時會發生頻散和衰減現象.研究地震波衰減和速度頻散特征對巖性分類、油氣識別，以及油氣田開發過程中的油氣藏管理都具有重要的價值(何潤發等，2020；凌云，2021；趙平起等，2021).由于地下介質條件復雜，地震波在地下含流體孔隙介質中傳播時的衰減和頻散特征一直是研究的熱點問題.

近幾十年來，國內外學者開展了大量巖石物理性質的研究，特別是孔隙介質流體的流動引起波的頻散衰減過程，發展出許多理論模型.

Biot(1955)建立了完全飽和孔隙介質中波傳播的動力學方程.在Biot理論中，黏性的孔隙流體和固體骨架之間的相互作用導致波的速度頻散和衰減.Geertsma和Smit(1961)應用Biot理論推導出縱波的速度頻散和衰減的近似解.然而，越來越多的研究發現，Biot理論的這種宏觀的流體流動機制很難解釋實際生產中地震波的高頻散和高衰減現象(Carcione et al.，2010).

為了突破Biot理論的局限性，學者們著眼于研究微觀尺度的流體流動機制.Mavko和Nur(1979)建立了噴射流模型來描述彈性波在流體部分飽和的裂縫/孔隙介質中的傳播，而且認為噴射流動會造成更強的衰減.歐陽芳等(2021)對經典噴射流模型進行了擴展，結合等效介質理論和孔隙結構模型，通過數值模擬研究了微觀孔隙結構下的速度頻散和衰減特性.

Dvorkin和Nur(1993)認為，在含有流體的巖石中，波的速度頻散和衰減除了受Biot機制的影響，還受到噴射流動機制的影響.他們將上述兩種機制統一起來，提出了BISQ模型，并推導了縱波相速度和品質因子的表達式，研究了滲透率等對縱波速度頻散和衰減的影響.

在Biot理論的框架下，White(1975a，b)首先引入中觀尺度的耗散機制，提出非均勻部分飽和的孔隙介質模型，也稱為斑狀飽和模型.Dutta和Odé(1979a，b)通過嚴格求解Biot方程組，研究了在鹽水充填的巖石中，衰減系數隨頻率、氣體飽和度等的變化.Carcione等(2003)深入研究了White模型的物性及流體變化對衰減和頻散的影響.王峣鈞等(2014)應用斑塊飽和巖石物理模型，從理論上分析了不同固結程度巖石中含氣飽和度對速度和衰減的影響，結果表明局部含氣儲層在地震頻帶內會發生頻散和振幅頻變效應.在此基礎上，他們建立了巖石物理量板，用于估計儲層的含氣飽和度.

Sun(2021)提出了一種BIPS(Biot-Patchy-Squirt)模型，用于描述非混相流體飽和裂縫孔隙彈性介質中的波的頻散衰減特性.研究結果表明，BIPS模型與實驗室數據吻合較好，這一發現將推動上述三種頻散衰減機制在波速預測中的潛在應用.

隨著研究的深入，學者們開始建立和研究多重孔隙介質模型.巴晶等(2012)建立了雙重孔隙結構模型，利用描述非飽和雙重孔隙介質地震波傳播的Biot-Rayleigh方程求解平面波解，得到縱波、橫波的相速度及逆品質因子的表達式，研究了非飽和巖石中的縱波頻散與衰減特征.通過對三個地區砂巖儲層的分析，發現在地震頻帶內縱波對儲層中含氣較為敏感，但對含氣飽和度的定量表征效率不高，低孔隙度砂巖的縱波頻散和衰減在地震頻帶更為顯著.此后，郭夢秋等(2018)進一步研究了基于此模型的含流體致密砂巖的縱波頻散和衰減特征.Ba等(2017)通過分析流體非均質性與巖石組構非均質性的疊加效應，提出了一種雙重雙孔隙理論.該模型描述了波在不同尺度下的具有組構非均質性的斑狀飽和巖石中的傳播.建模和實驗數據表明，巖石組構的非均質性尺度會造成波的多頻頻散和衰減效應，特征頻率的大小也與巖石組構的非均質性尺度有關.Zhang等(2020，2021)考慮到巖石組構具有多尺度的非均質性和自相似的分形特征.在此基礎上，他們建立了多孔介質模型，通過數值模擬研究了波的頻散和衰減特性與分形維數的關系.

在含油氣復雜儲層勘探時，不可避免的要遇到含裂縫和孔隙的儲層巖石，這些巖石孔隙結構復雜、孔隙流體性質多變，導致波場信號出現頻散衰減等特征.建立合理的裂縫-孔隙介質模型，是地震正演等技術的重要基礎.

Chapman等(2002)提出了在微觀尺度上建立孔隙彈性介質模型的方法，將模型設置成規則的立方體網格，每個網絡節點代表孔隙/裂縫空間.但是此模型并未建立孔隙介質滲透率與裂縫參數等之間的幾何關系.Tang(2011)、Tang等(2012)提出了含裂縫孔隙介質的波動方程模型，在推導過程中直接使用了Johnson等(1987)提出的動態滲透率模型，并沒有在滲透率與裂縫參數等影響因素之間建立關系，也沒有給出波頻散衰減與滲透率之間的關系.Song等(2020)假定裂縫的形狀可以是矩形的，建立了一個介質模型來估計飽和多孔巖石中包含一個隨機分布的、無限小厚度的矩形裂縫時的縱波衰減和速度頻散特征.

在天然巖石中，流體流動發生的空間是全局性的裂縫-孔隙網絡，而不是局部的單一孔隙或孔管束.Xiong等(2020)、熊繁升等(2021)基于巖石內部的裂縫連通網絡，用橢圓截面縱橫比的變化模擬從扁裂縫、軟孔隙到硬孔隙的多種情況，提出了具有橢圓形截面的三維裂縫/軟孔隙網絡模型，建立了宏觀可測觀量與裂縫參數之間的關系，并給出了滲透率的計算方法.該模型可以實現孔隙介質中跨尺度的整體建模，能更好地反映孔隙系統的整體效應.

考慮到熊繁升等(2021)建立的流體全飽和的三維裂縫/軟孔隙網絡模型沒有考慮孔隙的存在，因此本文在此模型的基礎上建立同時包含裂縫和孔隙的三維裂縫-孔隙網絡模型，并推導出滲透率的計算方法.基于三維裂縫/軟孔隙網絡模型和三維裂縫-孔隙網絡模型，運用體積平均的方法推導出基于滲透率的波動方程，利用平面波分析方法得到縱波頻散和衰減曲線的表達式，詳細研究了總孔隙度、裂縫孔隙度、裂縫縱橫比、裂縫數密度以及孔隙流體黏度對快縱波頻散衰減特征的影響，以及分析了慢縱波的頻散衰減特征.

1 三維孔隙網絡模型

1.1 三維裂縫/軟孔隙網絡模型

熊繁升等(2021)提出三維裂縫/軟孔隙網絡模型(圖1).模型中單個管道橫截面是縱橫比可變的橢圓形，可以模擬從圓形孔隙到狹窄裂縫的不同類型的裂縫/軟孔隙網絡空間.

圖1 三維裂縫/軟孔隙網絡模型Fig.1 3D fracture/soft pore network model

設單元立方體邊長為l，沿X、Y、Z方向上分別排布有M、N、L個橢圓截面管道，用R1和R2表示裂縫橢圓截面的長短軸半徑.定義如下參數：

(1)微裂縫的縱橫比：a=R2/R1.

(2)裂縫的孔隙度：

(1)

(3)單元體內的裂縫數密度：

(2)

熊繁升等(2021)計算三維裂縫/軟孔隙介質滲透率的步驟如下：

第一步：采用橢圓柱坐標系，基于質量守恒和NS方程推導出含不可壓縮牛頓流體在橢圓截面微管中的流量表達式.

第二步：利用微管兩端每個節點的流量守恒條件，得到全部網絡節點滿足的流量-壓力線性方程組.該方程組中未知量為各節點的壓力，入口和出口端的壓力邊界條件是方程的非齊次項.

第三步：求解上述流量-壓力線性方程組，得到各節點壓力，進而計算得到整個三維裂縫/軟孔隙網絡的流量.

第四步：結合達西定律計算得到巖石樣本的滲透率.

1.2 三維裂縫-孔隙網絡模型

在三維裂縫/軟孔隙網絡模型的基礎上，本文建立了同時包含孔隙和裂縫的三維裂縫-孔隙網絡模型(圖2).

圖2 三維裂縫-孔隙網絡模型Fig.2 3D fracture-pore network model

同樣地，設單元立方體邊長為l，沿X、Y、Z方向上分別排布有M、N、L個橢圓截面管道.節點上的球形表示孔隙，連接節點的橢圓柱體通道表示裂縫.用rpore表示孔隙半徑，Rasp表示平均孔隙半徑與平均裂縫截面長軸半徑之比.定義連接孔隙i、j的裂縫截面長軸半徑為

(3)

定義如下參數：

(1)微裂縫的縱橫比a：裂縫橢圓截面的長短軸半徑之比的倒數；

(2)總孔隙度：φ=φf+φpore；

其中，裂縫所占孔隙度為

(4)

孔隙空間所占孔隙度為

(5)

(3)單元體內的裂縫數密度：表達式同公式(2).

對于三維裂縫-孔隙網絡模型的滲透率計算，仍采用1.1節的推導方法.

對于同時含有裂縫和孔隙網絡的介質，考慮任意一段微管，微管兩端的網絡節點為i和j，節點處壓力分別是pi和pj，微管中壓力與流量的關系式為

Qij=cijpi+dijpj，

(6)

其中：

(7)

(8)

這里：

(9)

(10)

(11)

接下來，對各節點列流量平衡方程式，得到全部網絡節點滿足的流量-壓力方程組.求解該方程組得到各節點壓力，進而計算得到整個三維裂縫-孔隙網絡的流量.最后，結合達西定律計算得到巖石樣本的滲透率.

2 基于體積平均法的波動方程推導

2.1 體積平均法

體積平均法是推導波動方程的重要方法.Whitaker(1999)詳細介紹了采用體積平均法推導波動方程的原理和步驟.由于該方法通過對微觀物理量在單元體上進行積分平均得到宏觀物理量，因此可以考慮不同流體類型和孔隙空間的具體幾何特征，可以認為是一種更一般化的建模方法.

記單元體Ω的體積為V，孔隙空間由固體骨架的體積和孔隙介質中流體的體積構成，即：

V=Vs+Vf，

(12)

記ψf為孔隙介質中與流體有關的物理量，規定在流體的外部，ψf的值為0.定義在整個區域Ω上對ψf作體積平均的方法為

(13)

定義本征體積平均量：

(14)

根據以上定義式，可知：

(15)

2.2 波動方程推導

在下述的推導過程中：假設流體為牛頓流體，不考慮流體發生正應變過程中黏性的作用，僅考慮剪切黏性；考慮孔隙介質含一種固體骨架、一種流體(或兩種流體等效成一種).

用u、U分別表示固體、流體質點的位移向量；各物理量上方的一點表示對時間求導；以上下標中的s、f分別表示固體骨架、孔隙介質中所含的流體.

下面運用體積平均法來推導基于滲透率的三維孔隙網絡模型(三維裂縫/軟孔隙網絡模型與三維裂縫-孔隙網絡模型)的波動方程.

第一步，建立固體和流體的微觀運動方程(動量守恒方程).

固體骨架部分的微觀運動方程為

(16)

流體的微觀運動方程為

(17)

式中，ρs為固體骨架的質量密度；σ為固體應力；ρf，pf為等效后流體的密度、壓力；?表示張量積；s為流體應力.

第二步，給出固體部分和流體的本構方程.

固體部分的微觀本構方程為

(18)

式中，Ks為固體顆粒的體積模量；G為固體骨架的剪切模量；e為固體的體應變.

流體本構方程為：

(19)

式中，η為等效后流體的黏度.

第三步，對微觀方程做體積平均得到宏觀大尺度方程.

對流體微觀運動方程取體積平均后得到：

(20)

對固體微觀運動方程取體積平均后得到：

(21)

式中，κ為孔隙介質的滲透率，可分別根據1.1節和1.2節的計算方法確定.

第四步，推導出流固相互作用的流體和固體的本構關系.

固體和流體的本構關系為

(22)

式中，彈性常數aij可根據固體和流體的孔隙度、體積模量來進行計算.具體表示為

(23)

這里：

(24)

等效體積模量Kb的表達式為(Mavko and Nur，1978)：

(25)

式中，v是固體骨架礦物的泊松比；Nc為單位體內的裂縫數；R1i是第i條裂縫的橢圓截面長軸半徑長度；di是裂縫延伸長度；V是單元體的體積.

最后，將本構方程和體積平均后的動量守恒方程結合，得到最終形式的波動方程:

(26)

2.3 頻散/衰減表征

對波動方程兩端取散度，可得到縱波方程.對于各向同性孔隙介質，假設彈性波沿各方向傳播一致，利用平面波分析方法可得到縱波頻率-波數方程.將縱波頻率-波數方程的解記為v，根據如下定義，即可計算縱波頻散/衰減曲線：

(27)

3 頻散/衰減特征分析

波場在含流體的孔隙介質中傳播時會產生頻散和衰減現象，即波速隨著頻率發生改變，同時波的振幅也會發生變化.

下面基于三維裂縫/軟孔隙網絡模型和三維裂縫-孔隙網絡模型分別研究縱波頻散與衰減特征.

3.1 基于三維裂縫/軟孔隙網絡模型的頻散衰減特征

利用前面所述頻散衰減曲線的表達式，根據巖石物理參數可計算得到基于三維裂縫/軟孔隙網絡模型的頻率依賴的縱波速度和逆品質因子曲線.數值模擬所選取的模型、礦物和流體參數具體見表1.

表1 模型、礦物及流體參數表Table 1 Parameters of model, mineral and pore fluid

3.1.1 裂縫孔隙度的影響

圖3、圖4給出了不同裂縫孔隙度條件下，快縱波的頻散和衰減曲線.本算例中除裂縫孔隙度外，其余參數見表1.

圖3、圖4中四條曲線的裂縫孔隙度分別為0.001、0.005、0.02和0.06，對應的滲透率分別為7.00×10-17m2、1.57×10-16m2、3.14×10-16m2和5.44×10-16m2.從圖3可以看出，隨著裂縫孔隙度的減小，快縱波速度增大，起始速度由3116 m·s-1增大到4863 m·s-1.而且隨著裂縫孔隙度的增大，快縱波速度表現出更明顯的頻散現象.從圖4可以看出，隨著裂縫孔隙度的減小，逆品質因子曲線的峰值下降，特征頻率(逆品質因子曲線的峰值對應的頻率)向低頻方向移動.

圖3 快縱波頻散曲線隨裂縫孔隙度的變化Fig.3 Velocity dispersion changing with fracture porosity

圖4 快縱波衰減曲線隨裂縫孔隙度的變化Fig.4 Inverse quality factor versus fracture porosity

3.1.2 裂縫縱橫比的影響

圖5、圖6給出了不同裂縫縱橫比條件下，快縱波的頻散和衰減曲線.本算例中除裂縫縱橫比外，其余參數見表1.

圖5、圖6中四條曲線的裂縫縱橫比分別為0.02、0.04、0.08和0.16，對應的滲透率分別為5.44×10-16m2、3.07×10-15m2、1.73×10-14m2和9.60×10-14m2.通過計算可知，隨著裂縫縱橫比的下降，滲透率下降.這是一個合理的事實，因為隨著裂縫縱橫比的減小，裂縫趨于閉合.圖5顯示，隨著裂縫縱橫比的增大，快縱波速度增大.這主要是因為隨著縱橫比的增大，干巖石的可壓縮性變差，干巖石的體積模量Kb增大，從而引起快縱波速度的增大.圖6顯示，隨著裂縫縱橫比的增大，特征頻率向低頻方向移動.

圖5 快縱波頻散曲線隨裂縫縱橫比的變化Fig.5 Velocity dispersion versus aspect ratio of fracture

圖6 快縱波衰減曲線隨裂縫縱橫比的變化Fig.6 Attenuation changing with aspect ratio of fracture

3.1.3 裂縫數密度的影響

圖7、圖8給出了不同裂縫數密度條件下，快縱波的頻散和衰減曲線.本算例中除網絡節點數外，其余參數見表1.

圖7 快縱波頻散曲線隨裂縫數密度的變化Fig.7 Velocity dispersion varying with fracture density

圖8 快縱波衰減曲線隨裂縫數密度的變化Fig.8 Attenuation varying with fracture density

沿X、Y、Z方向上分別排布有3×3×3、4×4×4、5×5×5和6×6×6個孔隙，即裂縫數密度分別為54、144、300和540.計算得到對應的滲透率分別為1.40×10-16m2、3.11×10-16m2、5.44×10-16m2和8.39×10-16m2.通過計算可知，滲透率隨裂縫數密度的增加而增大.這表明隨著裂縫之間的連通程度升高，巖石滲透率逐步增大.圖7顯示，頻散幅值(速度極大值與極小值的差)對裂縫數密度的變化不敏感，但頻散曲線隨裂縫數密度的變化發生移動.在相同的頻率上，快縱波速度最大相差僅10 m·s-1.圖8顯示，隨著裂縫數密度的增大，特征頻率向低頻方向移動，而幾乎不影響逆品質因子曲線的峰值.

3.1.4 孔隙流體黏度的影響

下面研究孔隙流體的黏度對快縱波頻散衰減特征的影響.孔隙中的流體參數見表2，其余參數同表1.

表2 孔隙流體參數表Table 2 Parameters of pore fluid

從圖9、圖10可以看出，在此模型下的快縱波速度：油飽和砂巖明顯大于水飽和砂巖，水飽和砂巖明顯大于氣飽和砂巖.這主要是因為孔隙流體黏度增大時，孔隙流體與巖石之間的耦合作用加強，相當于提高了巖石的總體剛度，從而提高了快縱波的速度.但油飽和時的速度變化率相對于水飽和時變小.相較于孔隙充填其他流體，氣飽和時的快縱波頻散現象最不明顯，且逆品質因子曲線的峰值最小.

圖9 不同孔隙流體黏度下的快縱波頻散曲線Fig.9 Velocity dispersion curve as function of fluid viscosity

圖10 不同孔隙流體黏度下的快縱波衰減曲線Fig.10 Attenuation curve under different fluid viscosity

模型內不僅存在快縱波，還存在慢縱波.將表1中的巖石物理參數代入公式(27)，得到慢縱波的頻散衰減曲線，如圖11、圖12所示.

圖11 慢縱波頻散曲線Fig.11 Velocity dispersion curve of slow P-waves

圖12 慢縱波衰減曲線Fig.12 Attenuation curve of slow P-waves

3.1.5 慢縱波頻散衰減曲線

從圖11可以看出，在低頻范圍內(0～102Hz)慢縱波速度基本為0，且在0～108Hz頻率范圍內發生頻散，速度隨頻率的增加而增大，在高頻范圍內穩定在大約740 m·s-1.從圖12可以看出，慢縱波逆品質因子曲線在0～106Hz的頻率范圍內為水平線，后隨著頻率的增加而減小，呈直線下降趨勢.

3.2 基于三維裂縫-孔隙網絡模型的頻散衰減特征

根據式(3)知，裂縫橢圓截面長軸半徑與裂縫兩端的孔隙半徑有關.因此，裂縫的孔隙度φf和孔隙空間的孔隙度φpore不是互相獨立的.

在實際的計算中，給定φpore，首先根據式(5)計算得到孔隙半徑，后根據式(3)計算得到裂縫截面的長軸半徑，再根據裂縫縱橫比得到裂縫的短軸半徑.將上述計算得到的裂縫截面的長軸半徑代入式(4)，計算得到φf.

為了滿足預設的總孔隙度φ，可能需要對孔隙半徑和裂縫半徑進行重新處理.因此，在三維裂縫-孔隙網絡模型中，總孔隙度φ與預設結果相同，但是φpore與φf可能與預設結果不同.

3.2.1 總孔隙度的影響

除裂縫孔隙度外，仍采用表1中的巖石物理參數.圖13、圖14給出了總孔隙度分別為0.001、0.005、0.02和0.06時的快縱波頻散與衰減曲線.通過計算得到對應的滲透率為3.12×10-17m2、7.00×10-17m2，1.40×10-16m2和2.22×10-16m2.從圖中可以看出，隨著總孔隙度的減小，快縱波速度增大，起始速度由3559 m·s-1增大到4367 m·s-1，逆品質因子曲線的峰值下降，特征頻率向低頻方向移動.

圖13 快縱波頻散曲線隨總孔隙度的變化Fig.13 Velocity dispersion curve as function of total porosity

圖14 快縱波衰減曲線隨總孔隙度的變化Fig.14 Attenuation curve as function of total porosity

3.2.2 裂縫參數的影響

從建模過程中可以發現，三維裂縫/軟孔隙網絡模型可以認為是三維裂縫-孔隙網絡模型的特例，即當孔隙空間所占孔隙度為0時，三維裂縫-孔隙網絡模型退化為三維裂縫/軟孔隙網絡模型.若固定孔隙空間所占孔隙度為0，分析裂縫參數對快縱波頻散衰減特征的影響，即為文章3.1節的內容.

3.2.3 慢縱波頻散衰減曲線

將表1中的巖石物理參數代入公式(27)，得到慢縱波的頻散衰減曲線.

從圖15、圖16可以看出，在該模型下慢縱波頻散與衰減的趨勢與三維裂縫/軟孔隙網絡模型相似，但頻散的幅值更大，發生頻散的頻率范圍更寬，在高頻范圍內慢縱波的速度穩定在大約1214 m·s-1.

圖15 慢縱波頻散曲線Fig.15 Velocity dispersion curve of slow P-waves

圖16 慢縱波衰減曲線Fig.16 Attenuation curve of slow P-waves

4 結論

本文基于三維裂縫/軟孔隙網絡模型和三維裂縫-孔隙網絡模型，詳細研究了總孔隙度、裂縫孔隙度、裂縫縱橫比、裂縫數密度、流體黏度對快縱波頻散衰減特征的影響，以及分析了慢縱波的頻散衰減特征.

基于三維裂縫/軟孔隙網絡模型，可以發現：

(1)隨著裂縫孔隙度的增大，快縱波速度表現出更明顯的頻散現象.隨著裂縫孔隙度的減小，快縱波速度增大，逆品質因子曲線的峰值下降，特征頻率(逆品質因子曲線的峰值對應的頻率)向低頻方向移動.

(2)隨著裂縫縱橫比的增大，快縱波速度增大，特征頻率向低頻方向移動，而頻散幅值(速度極大值與極小值的差)和逆品質因子曲線的峰值的大小未表現出規律性的變化.

(3)裂縫數密度的變化幾乎不會影響快縱波速度、頻散幅值和逆品質因子曲線的峰值，但對特征頻率影響顯著.隨裂縫數密度的增大，特征頻率向低頻方向移動.

(4)當孔隙充填油、水、氣時，油飽和砂巖的快縱波速度明顯大于水飽和砂巖，水飽和砂巖明顯大于氣飽和砂巖.氣飽和時的快縱波頻散現象最不明顯，且逆品質因子曲線的峰值最小.

基于三維裂縫-孔隙網絡模型，可以發現：

在此模型下，總孔隙度、裂縫參數等對快縱波頻散衰減特征的影響以及慢縱波的頻散衰減趨勢與三維裂縫/軟孔隙網絡模型相似.

總的來說，孔隙度的變化主要影響逆品質因子曲線的峰值的大小；裂縫數密度主要控制速度顯著變化的范圍；裂縫縱橫比的變化對縱波速度和特征頻率都影響顯著.