數字巖心寬頻帶動態應力應變模擬方法及其對含裂隙致密巖石頻散和衰減特征的表征

朱偉， 趙巒嘯， 王一戎

1 長江大學油氣資源與勘探技術教育部重點實驗室， 武漢 430100 2 長江大學地球物理與石油資源學院， 武漢 430100 3 同濟大學海洋地質國家重點實驗室， 上海 200092 4 同濟大學海洋與地球科學學院， 上海 200092

0 引言

非均質性是地下多孔巖石的最基本表征之一.在漫長的地質年代，由于沉積和成巖作用的差異，以及構造運動和地質過程中的隨機擾動，巖石往往在各個尺度上表現出不同程度的非均質性.這些非均質性表現為巖石的物質組成、孔隙類型、流體分布、裂隙發育、輸運性質等空間分布的不均勻性(Zhao et al.，2020).當彈性波在含流體非均勻孔隙介質中傳播時，地層的非均質性不僅影響彈性波的傳播速度，而且由于波致孔隙壓力差異造成的流體流動還會引起彈性波的頻散和衰減效應(Müller et al.，2010).孔隙壓力弛豫在不同類型的巖石中造成的彈性波頻散和衰減特征已經被大量的實驗數據和不同頻段的地球物理數據所證實(未睍等，2015；Subramaniyan et al.，2015；Chapman et al.，2019；Ma et al.，2018；Borgomano et al.，2019；李闖等，2020；Li et al.，2020).研究彈性波在不同頻段的速度和衰減特征一方面可以指導綜合利用多尺度的地球物理數據進行地層非均質性的刻畫，另外一方面為利用地震耗散相關的屬性進行流體和裂隙儲層等預測提供物理基礎.同時，這些研究對于非均質儲層的地震表征、四維地震的流體監測、二氧化碳封存、地熱資源的勘探開發都具有重要的價值(Sams et al.，1997；Zhao et al.，2015；Chen and Innanen，2018;Chen，2020).

近幾十年來，很多學者提出大量理論模型來刻畫波致流體流動引起的頻散和衰減效應(Biot，1962；White，1975；Dvorkin et al.，1995；Gurevich et al.，1997；Chapman et al.，2002；Pride et al.，2004；Ba et al.，2008, 2017；Tang，2011；巴晶等，2012；Yao et al.，2015；Zhao et al.，2017；Papageorgiou and Chapman，2017；Jin et al.，2018；趙正陽等，2019)，這些理論大多基于對巖石非均質特征(比如孔隙結構或者分布形式)的一些簡化和假設，聚焦于探討不同的地質和物理參數對頻散和衰減特征的影響.相對而言，數字巖石物理技術則可以很好地彌補這種不足，可以模擬各種接近實際復雜非均質多孔巖石的孔隙結構、裂隙形態、流體分布，從而可以全面準確刻畫復雜非均質多孔巖石中波致流體流動引起的頻散和衰減效應.這些數字巖心的尺度雖然比實際地層要小很多，但其揭示的基本物理過程和規律卻對綜合利用多頻段地球物理數據進行多尺度地質特征刻畫有重要的啟示意義.

數字巖心是數字巖石物理的模型基礎，可以用數學重建方法或物理實驗方法建立.在一定空間尺度上，數學重建方法建立的數字巖心與真實巖心的孔隙結構具有統計相似性，物理實驗方法建立的數字巖心與真實巖心的孔隙結構十分接近(朱偉和單蕊，2014).在數字巖心中，固體骨架和孔隙流體是相互獨立的介質，通過邊界連接在一起.質點的介質類型是固體和流體中的一種，而不是多相孔隙介質模型中固體和流體共存的情況.數字巖石物理中的彈性/黏彈性模擬技術主要包括四種：1)靜力學模擬(Arns et al.，2002；Zhang et al.，2016)；2)透射波模擬(Saenger and Shapiro，2002；Saenger et al.，2011；印興耀等，2016；Li et al.，2019)；3)準靜態應力應變模擬(Quintal et al.，2016，2019；Alkhimenkov et al.，2020a，2020b；Lissa et al.，2020)；4)動態應力應變模擬(Zhang and Toks?z，2012；Zhu et al.，2017；Das et al.，2019；朱偉等，2020).Zhang和Toks?z(2012)首次提出了數字巖心動態應力應變模擬方法，計算了砂巖數字巖心在頻率為100 kHz時的速度和衰減.Zhu等(2017)建立了一種沿不同方向計算數字巖心的速度和衰減的動態應力應變模擬方法.Das等(2019)利用動態應力應變模擬技術分析頻散和裂隙扁率、頻散與流體黏滯系數的關系.朱偉等(2020)利用動態應力應變模擬方法研究砂巖CT數字巖心內部的位移場和孔隙壓力的空間和時間變化.

現有的數字巖心動態應力應變模擬為單頻率模擬方法，通過使數字巖心發生簡諧型應變計算某一頻率的速度和衰減.單頻率模擬可以提取任意模擬時刻數字巖心的流體位移和孔壓變化，并與速度、衰減進行關聯分析，但須在不同頻率進行多次模擬才能得到頻散和衰減曲線的基本形態.為了不漏掉重要的細節特征，模擬頻率必須仔細選擇，增加了數值模擬的復雜性.

因而，為了較為全面準確刻畫非均質多孔巖石在全頻段(地震—超聲)的頻散和衰減特征，需要發展相應的數字巖心模擬技術.寬頻帶動態應力應變模擬期望給數字巖心加載一個恒定應變，模擬數字巖心內部應力松弛過程，通過一次模擬計算一個頻帶內的頻散和衰減曲線.對地震巖石物理學而言，寬頻帶的頻率范圍是零至某個最大頻率，模擬時只需考慮最高截止頻率即可.在寬頻帶模擬中也可以提取任意時刻的流體位移和孔壓變化，但不能與某個頻率簡單聯系起來.

裂隙在多孔巖石尤其是致密地層中普遍存在.雖然裂隙在巖石中只占很小的體積百分比，但對巖石的彈性和滲透性有極大的影響.在一些致密儲層中，從地震上找到裂隙發育區的“甜點”成為致密儲層開發的關鍵，理解含裂隙致密巖石的地震彈性和衰減特征是裂隙屬性地震評價的基礎.當地震波通過含裂隙的致密巖石時，裂隙由于可壓縮性要大于背景孔隙(如粒間孔、溶蝕孔洞等)，因而會引起地質流體在裂隙和背景孔隙之間交流，并引起彈性波的頻散和衰減，一般被稱為擠噴流效應.

本文首先詳細介紹寬頻帶動態應力應變模擬的原理；然后將該方法用于刻畫含裂隙致密巖石擠噴流效應引起的速度頻散和衰減特征，并通過數字巖心模擬較為系統的揭示了致密地層中控制擠噴流效應的主控物理因素，也證明了該方法在表征復雜非均質多孔巖石波致流體流動引起的頻散和衰減特征的適用性.

1 基本原理

1.1 邊界條件

當前地震勘探仍以縱波勘探為主，對數字巖心設置如圖1所示的邊界條件，可計算縱波沿x軸方向傳播時的頻散和衰減.在數字巖心的上、下邊界應用周期性邊界條件，相當于將數字巖心的上、下兩個邊界連接在一起，即如果波從上邊界出射，則波將從下邊界上相同的水平位置以相同的方向入射，反之亦然.由于數字巖心在垂直方向受限，這使計算水平方向的縱波傳播性質變得容易.

圖1 數字巖心的邊界條件示意圖Fig.1 Setup of boundary condition on digital rock

數字巖心的左邊界靜止不動，右邊界受到外部動態作用，數字巖心發生水平方向的受迫運動.在t=0時刻，數字巖心中質點的速度和位移均為零，即

vx(x,z,t=0)=0,

vz(x,z,t=0)=0,ux(x,z,t=0)=0,

uz(x,z,t=0)=0.

(1)

式中，vx和vz是質點在x和z方向的速度分量，ux和uz是質點在x和z方向的位移分量.

當t≥0時，數字巖心左邊界質點的速度分量為

vx(0,z,t)=0,

vz(0,z,t)=0.

(2)

當t≥0時，數字巖心右邊界質點的速度分量為

vx(L,z,t)=A·f(t),

vz(L,z,t)=0.

(3)

式中，L表示數字巖心的邊長；f(t)是t的函數，f(0)=0；A是調節應變大小的常數，使數字巖心在x方向滿足小應變假設.

函數f(t)的形態是寬頻帶模擬和單頻率模擬的核心區別.在單頻率模擬(Zhu et al.，2017)中，f(t)是某個頻率的正弦函數.在寬頻帶模擬中，f(t)相當于一個低通信號，在通帶內振幅近似為常數，在壓制帶內振幅近似為零，模擬得到通帶內的頻散和衰減曲線.函數f(t)的一個可能實現如圖2a所示，類似于一個延遲尖脈沖函數，它表示數字巖心右邊界質點的速度分量vx隨時間的變化，而右邊界質點的位移分量ux隨時間的變化可以用圖2b表示.質點的位移是速度的時間積分，因此對圖2a中的曲線進行時間積分就得到圖2b中的曲線(圖2中曲線幅度為歸一化值).當模擬時間增大時，數字巖心右邊界質點的位移分量ux趨于某一穩定值.因此，在右邊界上的外力作用下，數字巖心的宏觀應變逐漸趨于定值，而數字巖心內部的應力場逐漸平均化，發生應力松弛.圖2c是f(t)的振幅譜，具有低通性質，在通帶內振幅值基本穩定.若f(t)越接近尖脈沖函數，則其通帶越寬.

圖2 歸一化的f(t)、f(t)的積分和f(t)的振幅譜(a) 函數f(t)； (b) f(t)的積分; (c) f(t)的振幅譜.Fig.2 Normalized f(t), integral of f(t) and amplitude spectrum of f(t)(a) The function f(t); (b) The integral of f(t); (c) The amplitude spectrum of f(t).

1.2 波場模擬

假設數字巖心的骨架由各向同性的線彈性介質組成，其本構方程為

(4)

式中，λ和μ是拉梅常數；σxx，σzz和σzx表示應力；exx，ezz和ezx表示應變.

假設數字巖心飽和牛頓流體，牛頓流體的本構關系為

(5)

式中，ηλ和ημ表示黏滯系數，在模擬中令ηλ=ημ=η.

用位移和應力表示的運動方程為

(6)

式中，ρ表示密度.

求解方程組(4)、(5)和(6)的方法是旋轉交錯網格有限差分法(Saenger et al.，2000).圖3表示旋轉交錯網格的網格設置，位移和速度分量位于網格頂點，應變、應力、彈性模量、黏滯系數和密度定義在網格中心.

圖3 旋轉交錯網格示意圖Fig.3 Sketch map of rotated staggered grid

若采用2階精度的旋轉交錯網格有限差分格式，則速度和位移的有限差分可直接用下面的公式表示.

(7)

式中，ψi可表示網格頂點的位移和速度，下標i=1,2,3,4表示網格的4個頂點的標號.

公式(7)可視為兩個中心差分格式的平均，用網格的4個頂點的速度和位移計算它們在網格中心的有限差分值.應變和應變率分別是位移和速度的空間導數.故應變、應變率、彈性模量和黏滯系數都定義在網格中心，更新應力時不需要對彈性模量和黏滯系數插值，流固邊界附近的應力穩定.當更新網格頂點的速度和位移時，需要使用以該頂點為公共頂點的4個網格來計算應力的有限差分.由于密度定義在網格中心，網格頂點的密度需要由周圍4個網格的密度進行平均.密度插值對迭代算法是穩定的.因此，在旋轉交錯網格有限差分中不需要對流固邊界做特別處理.

1.3 縱波的速度與逆品質因子

在模擬過程中計算數字巖心的平均正應變和正應力，在模擬結束后計算應變率和應力率.數字巖心的復縱波模量是應力率和應變率的傅立葉變換的比值.數字巖心的縱波速度和逆品質因子由復縱波模量和密度計算.上述過程可以用公式(8)表示.

(8)

1.4 精度驗證

本文采用均勻牛頓流體模型驗證模擬方法的精度.牛頓流體的體積模量為2.3 GPa， 黏滯系數為1.0 Pa·s，密度為1000 kg·m-3.模型網格的長度為1.0×10-6m，模擬時間步長為1.0×10-10s，模擬時間長度為1.0×10-3s.數值模擬計算的頻散和逆品質因子如圖4a和圖4c所示.由公式(5)可知，牛頓流體的復縱波模量為M(ω)=λ+i3ωη，可利用公式(8)的后三式計算頻散和逆品質因子的解析解，顯示于圖4a和圖4c中.

由圖4可知，1)頻散和逆品質因子的數值解與解析解基本重合，誤差非常小(圖4b和圖4d)；2)對于均勻牛頓流體而言，速度和逆品質因子隨頻率單調增加，但在圖4的頻率范圍內頻散和衰減程度均很小.

圖4 牛頓流體的頻散和逆品質因子(a) 速度頻散曲線； (b) 速度的數值解與解析解之差； (c) 逆品質因子曲線； (d) 逆品質因子的數值解與解析解之差.Fig.4 Velocity dispersion and inverse quality factor of Newtonian fluid(a) Velocity dispersion curve; (b) Velocity difference between simulation and theory; (c) Inverse quality factor curve; (d) Inverse quality factor difference between simulation and theory.

2 含裂隙致密巖石的頻散和衰減特征數值模擬

這一部分將數字巖心寬頻帶動態應力應變模擬方法用于表征含裂隙致密巖石的擠噴流效應，一方面旨在細究控制含裂隙致密地層擠噴流效應的主控物理因素，另一方面也證明該方法在刻畫復雜非均質多孔巖石由于波致流體流動引起的速度頻散和衰減特征的有效性.當巖石受到擠壓或拉伸時，裂隙中的流體壓力變化較大，在孔隙中形成較大的壓力梯度，驅動孔隙壓力在“較軟”的裂隙和較硬的孔隙之間擴散，導致波的頻散和衰減.針對含裂隙致密巖石的物理特征，本節建立了4個二維數字巖心，用寬頻帶動態應力應變數值模擬方法計算數字巖心的頻散和逆品質因子曲線，并分析其特征.

2.1 數字巖心構建

為了刻畫擠噴流效應，數字巖心包含若干相互獨立的裂隙—圓孔系統.裂隙為直裂隙，走向垂直外部擠壓(拉伸)方向，一端與圓孔相連，另一端封閉，見圖5.數字巖心a和b的邊長為1.0×10-3m，孔隙度分別約為0.020和0.051.裂隙的長度為8.8×10-5m，寬度為3×10-6m，扁率(也稱縱橫比，定義為矩形裂隙的寬度與長度的比值)約為0.034，圓孔半徑為1.5×10-5m.數字巖心c和d的寬度為0.5×10-3m，高度為1.0×10-3m，孔隙度分別約為0.024和0.052.裂隙的長度為8.8×10-5m，寬度為1.5×10-6m，扁率約為0.017，圓孔半徑為1.5×10-5m.

這里采用硬幣狀裂隙的近似公式Cr=3φc/4πα計算裂隙密度，φc為裂隙孔隙度，α為裂隙扁率(Hudson，1981；Mavko et al.，2009).數字巖心a、b、c和d的裂隙密度分別約為0.035，0.089，0.048和0.104.

圖5 數字巖心數字巖心a和b的邊長為1.0×10-3 m，孔隙度分別為0.020和0.051.裂隙長度約8.8×10-5 m，裂隙寬度約3×10-6 m，扁率約0.034，圓孔半徑約為1.5×10-5 m.數字巖心c和d的寬度為0.5×10-3 m，高度為1.0×10-3 m，孔隙度分別為0.024和0.052.裂隙長度約8.8×10-5 m，裂隙寬度約1.5×10-6 m，扁率約為0.017，圓孔半徑約為1.5×10-5m.Fig.5 Digital rocks Side length of digital rock a and b is 1.0×10-3 m, their porosity are 0.02 and 0.051. The length and width of cracks are 8.8×10-5 m and 3×10-6 m, the aspect ratio is 0.034. Radius of pores is 1.5×10-5 m. The width and height of digital rock c and d are 0.5×10-3 m and 1.0×10-3 m, their porosity are 0.024 and 0.052. The length and width of cracks are 8.8×10-5 m and 1.5×10-6 m, the aspect ratio is 0.017, Radius of pores is 1.5×10-5 m.

這4個數字巖心的骨架由同一固體介質構成，孔隙被同一流體飽和.固體與流體的參數見表1.綜合考慮模擬效果、裂隙扁率、模型大小和計算負擔，我們對流體設計了較大的黏滯系數(朱偉等，2020).

表1 數字巖心骨架固體和孔隙流體的黏彈性參數表Table 1 Viscoelasticity of solid matrix and pore fluid in digital rock

數值模擬采用基于正方形網格的有限差分法，在對數字巖心離散化后，孔隙和骨架的界面局部呈現鋸齒狀，這可能會影響波致流體流動.但是，數字巖心的裂隙為直裂隙，裂隙面與網格邊界走向一致，離散后的裂隙表面仍然是光滑的，不存在鋸齒狀的情況.當數字巖心受到垂直裂隙面方向的擠壓(拉伸)時，強烈的流體運動發生在裂隙中.所以，裂隙中的擠噴流基本不受鋸齒狀孔隙邊界的影響.

2.2 模擬測試

當裂隙飽和黏性流體時，需要對裂隙寬度進行充分采樣以便以適當的精度刻畫波致流效應.我們以數字巖心b為例進行網格剖分和數值模擬試驗，確定離散裂隙寬度所需的網格數量.首先，將裂隙寬度離散為3個網格；然后，以此為基礎，逐步加密網格(即減小網格大小)，分別得到裂隙寬度為6個和12個網格的模型.如此，我們得到3個離散的數字巖心網格模型，它們的大小和結構完全相同，只是網格大小和數量不同.裂隙寬度分別為3個、6個和12個網格.

在數值模擬中，3個模型采用相同的模擬時間步長和總時長，得到的頻散和逆品質因子曲線如圖6所示.函數f(t)的通帶范圍為0～2 MHz，當頻率高于1 MHz后，由于波場散射逐漸明顯，平均應力曲線中包含了較多擾動，頻散和逆品質因子曲線的高頻干擾比較嚴重，故在圖6中頻率上限為1 MHz.

由圖6可知，這3個模型的頻散和逆品質因子曲線的形態基本一致，但特征頻率發生變化.當裂隙寬度為12個和6個網格時，特征頻率非常接近.當裂隙寬度為3個網格時，特征頻率向低頻方向略有移動.

圖6 數字巖心b的頻散和逆品質因子曲線曲線3-grids、6-grids和12-grids分別對應裂隙寬度用3個、6個和12個網格離散.Fig.6 Dispersion and attenuation of digital rock b Curves named with 3-grids, 6-grids and 12-grids corresponding to crack width sampled by 3, 6 and 12 grids, respectively.

離散模型的網格越小，模擬結果的精度越高，這是數值模擬的一個基本特征.但是，網格減小會使網格數量增加、模擬步長減小，從而大幅提高計算工作量.因此，在精度合理的前提下，網格大一點對完成數值模擬有利.在本文中，我們認為用6個網格對裂隙寬度進行離散比較合適.進一步加密網格可以提高精度，但效果逐漸不明顯.

與網格大小有關的另一不利影響因素是數值頻散.但當波長與網格邊長的比值較大時，數值頻散可以勿略.當數字巖心的裂隙寬度分別為3個，6個和12個網格時，網格邊長分別為1.0×10-6m，0.5×10-6m和0.25×10-6m.彈性波在骨架礦物和孔隙流體中的縱、橫波速度分別約為6008 m·s-1，4075 m·s-1和1517 m·s-1.以最高頻率1 MHz計算的最短波長約為1.517×10-3m.波長與網格邊長的最小比值約為1517，即最短波長遠遠大于最大網格邊長，在此條件下數值頻散非常小，不影響擠噴流的頻散效應.

2.3 模擬結果

根據2.2小節的試驗結果，我們采用6個網格對裂隙寬度進行離散.相應地，數字巖心a和b的網格大小為0.5×10-6m，數字巖心c和d的網格大小為0.25×10-6m.兩組數字巖心的模擬時間步長分別為0.25×10-10s(a和b)和0.1×10-10s(c和d).數字巖心c和d的裂隙扁率比數字巖心a和b小.根據擠噴流理論推測數字巖心c和d中流體壓力平衡過程需要更長的時間.兩組數字巖心的模擬時間長度分別為1.25×10-4s(a和b)和5.0×10-4s(c和d).

圖7為上述4個數字巖心的模擬結果，圖7a和圖7b分別為頻散和逆品質因子.在圖7a、圖7b中，曲線a，b，c和d分別表示數字巖心a，b，c和d的模擬結果.

圖7 數字巖心的模擬結果(a) 速度; (b) 逆品質因子.曲線a、b、c、d分別表示數字巖心a、b、c、d的模擬結果，φa，φb，φc和φd表示總孔隙度，Cra，Crb，Crc和Crd表示裂隙密度.Fig.7 Simulation results of digital rocks(a) Velocity; (b) Inverse quality factor. Curve-a, -b, -c and -d are results of digital rock-a, -b, -c and -d respectively.φa,φb,φc and φd represent total porosity.Cra, Crb, Crc and Crd represent crack density.

頻散曲線表現的特征為：數字巖心a和c的速度明顯要比數字巖心b和d高，這主要是由于數字巖心b和d的裂隙密度較高造成的.裂隙是一種體積和扁率非常小，但對巖石彈性影響非常大的孔隙類型.這4個數字巖心中的裂隙是定向排列的，垂直裂隙面方向的波速受到裂隙的影響最大.當頻率小于103Hz時，數字巖心a與c的速度基本不隨頻率變化，這主要是由于在此頻段內，裂隙內流體近似處于松弛狀態；而數字巖心c比a的速度稍低，這主要還是由于數字巖心c的孔隙度比a稍高，所以數字巖心c在松弛階段的速度稍低.相較與數字巖心a，數字巖心c開始出現明顯頻散效應的頻率較低，這主要是由于數字巖心c的裂隙扁率較小，因此特征頻率較低.這符合一般的物理認識，也就是扁而長的裂隙往往需要更長的時間進行孔隙壓力弛豫，因此其衰減的特征頻率也會較低.而且由于頻散效應，數字巖心c的縱波速度很快超過了數字巖心a，這主要是因為數字巖心c稍高的裂隙密度造成的，這樣有更多的流體在較軟的裂隙和較硬的圓孔之間交流，引起了較大的頻散效應.類似于數字巖心c和a，數字巖心d的頻散效應明顯高于數字巖心b，并且特征頻率出現在較低的頻段.在低頻段數字巖心d的速度稍稍高于數字巖心b的速度，但數字巖心d的孔隙度卻略高于數字巖心b.這可以用孔隙結構的變化來解釋.

如圖7b所示，數字巖心的衰減特征與速度頻散有較好的對應關系.可以清楚地看到，數字巖心b與a、數字巖心d與c分別具有近似一樣的特征頻率，這也證明了在其他物理參數一樣的情況下，裂隙扁率對特征頻率的控制作用.而數字巖心b和d的衰減幅度明顯要高于數字巖心a和c，這也一定程度上反映了裂隙密度對衰減幅度的一階控制作用.另外，在特征頻率的低頻側，逆品質因子下降速率較快；而在特征頻率的高頻側，逆品質因子下降速率較慢，而且下降速率變化，下降速率快、慢交替出現.這個現象在圖7b(因為縱軸為逆品質因子的對數)中不易發現，但在圖6中比較明顯.

3 討論

3.1 單頻率與寬頻帶動態應力應變模擬的關系

單頻率模擬計算某一頻率的速度和逆品質因子，可提取任意模擬時刻孔隙流體壓力變化的快照進行分析.寬頻帶模擬可以計算一個頻帶范圍內的頻散和逆品質因子曲線.兩者聯合分析的基礎是它們的結果具有較好的一致性.圖8展示了數字巖心a和b的單頻率和寬頻帶動態模擬結果的對比圖，其中單頻率模擬的頻率為50 kHz、65 kHz和75 kHz.可以清楚地看到，在這三個頻率處，兩種方法的結果幾乎一致，說明兩種方法的一致性較好.

圖8 單頻率和寬頻帶模擬結果比較(a) 速度; (b) 逆品質因子.數字巖心為a和b，單頻率模擬的頻率為50 kHz、65 kHz和75 kHz.圖例中a-band和b-band表示數字巖心a和b的寬頻帶模擬結果，a-single和b-single表示數字a和b的單頻率模擬結果.Fig.8 Comparison of the results of single-frequency and broadband simulation (a) Velocity; (b) Inverse quality factor. The models are digital rock-a and -b, frequencies of the single-frequency simulation are 50 kHz, 65 kHz and 75 kHz. In legend, a-band and b-band represent the broadband simulation results of digital rock-a and -b, a-single and b-single represent the single-frequency simulation results of digital rock-a and -b.

3.2 含裂隙致密巖石的頻散和衰減特征的控制因素分析

含裂隙致密巖石的特征頻率與裂隙扁率的關系符合擠噴流模型的基本特征，即裂隙的扁率減小，則其對應的特征頻率向低頻方向移動.噴流模型的特征頻率(Saenger et al.，2011)為

(9)

式中，Km為基質體積模量，α為裂隙扁率，η為流體的黏滯系數.若特征頻率的變化僅由裂隙扁率引起，那么裂隙扁率分別為0.034和0.017的兩組數字巖心的特征頻率的比值為8.實際兩組數字巖心的特征頻率的平均比值約為7.數字巖心模擬的特征頻率比值與擠噴流理論模型計算的比值接近，驗證了裂隙扁率是擠噴流頻散效應特征頻率的主控因素，而兩者差異可能與裂隙-圓孔系統的物理描述、空間分布及其彈性相互作用都有關系.

圖9 衰減峰值與裂隙密度的關系Fig.9 Relationship between peak of inverse quality factor and crack density

圖9顯示了10個數字巖心(含數字巖心a、b、c和d)的逆品質因子的峰值與裂隙密度的關系.折線相連的實心圓表示5個裂隙扁率為0.017的數字巖心.折線相連的正方形表示5個裂隙扁率為0.034的數字巖心.由圖9可知，在圖示的裂隙密度范圍內，當裂隙扁率相同時，衰減峰值與裂隙密度總體上具有近似線性的正相關性，這與一些擠噴流理論模型的結論相似(Tang，2011；Yan et al.，2014)，但也受到孔隙結構的影響而產生局部變化.不同的裂隙扁率對衰減幅度-裂隙密度的關系有顯著的影響.裂隙密度相近條件下，狹長的裂隙(扁率較小)引起的衰減幅度更高.

從圖7可知，逆品質因子在低于特征頻率的一側隨著頻率的減小下降較快，而在高于特征頻率的一側隨著頻率的增加下降較慢，這應該主要是由于擠噴流效應造成的(Gurevich et al.，2010；Alkhimenkov et al.，2020a).另外，考慮到即使頻率為1 MHz時，波長與裂隙-圓孔系統總長度的比值可達40，根據頻散程度與波長-散射體尺度比值的關系，此時彈性散射引起的頻散和衰減程度是比較小的(Mavko et al.，2009).所以，本文動態應力應變數值模擬的結果基本上是擠噴流效應的反映.彈性散射對頻散和衰減的作用可能在更高的頻帶呈現，但前提是需要對平均應力曲線進行適當的壓噪處理.

3.3 對縫洞型儲層地震勘探的啟示

本例中頻散和衰減主要發生在聲波—超聲波頻帶內，也就是說對于含微觀尺度裂隙的巖石來說(<1 cm)，微裂隙與基質孔之間的擠噴流效應對測井數據的影響可能是顯著的.由于考慮到計算成本的問題，數字巖心的尺度較小，但所研究的模型對于大尺度的裂隙儲層來說也具有重要的啟示意義.比如塔里木盆地深層碳酸鹽巖儲層斷裂-溶蝕系統非常發育，因此中大尺度的裂隙往往與溶蝕洞構成復雜的儲集空間.而正如本研究所揭示的那樣，當地震波通過該類儲層時，中大尺度的裂隙與溶蝕洞之間同樣可以發生孔隙壓力的弛豫或者流體的交換，因此一樣會引起地震頻散和衰減.這里的數值模擬表明裂隙與孔洞間的流體交換所致的頻散和衰減的特征頻率與裂隙的扁率密切相關.也就是裂隙的張開度與裂隙長度的比值足夠低時，其特征頻率很可能發生在勘探地震的主頻內，從而對利用地震數據進行深層斷溶儲層的刻畫有重要啟示意義.

4 結論

為了刻畫復雜非均質多孔巖石在跨頻段(勘探地震、聲波測井、超聲波)的頻散和衰減特征，本文提出了基于數字巖心的寬頻帶動態應力應變模擬方法.該方法以固液耦合的波動方程為基礎，模擬數字巖心的應力松弛過程，計算數字巖心在目標頻帶范圍內的頻散和衰減特征.通過將該方法成功用于表征含裂隙致密巖石的擠噴流效應，證明了該方法可以用于有效刻畫復雜非均質多孔巖石波致流體流動引起的頻散和衰減特征.含裂隙-圓孔系統的數字巖心的寬頻帶動態應力應變數值模擬結果顯示，裂隙的扁率對速度頻散和衰減的特征頻率有很強的控制作用，較小的裂隙扁率對應著特征頻率向低頻移動；頻散和衰減的強度則與裂隙密度和裂隙扁率都有關系，較高的裂隙密度和較小的裂隙扁率對應著較強的衰減幅度.