朱偉， 趙巒嘯， 王一戎

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

0 引言

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

近幾十年來，很多學(xué)者提出大量理論模型來刻畫波致流體流動引起的頻散和衰減效應(yīng)(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)，這些理論大多基于對巖石非均質(zhì)特征(比如孔隙結(jié)構(gòu)或者分布形式)的一些簡化和假設(shè)，聚焦于探討不同的地質(zhì)和物理參數(shù)對頻散和衰減特征的影響.相對而言，數(shù)字巖石物理技術(shù)則可以很好地彌補這種不足，可以模擬各種接近實際復(fù)雜非均質(zhì)多孔巖石的孔隙結(jié)構(gòu)、裂隙形態(tài)、流體分布，從而可以全面準(zhǔn)確刻畫復(fù)雜非均質(zhì)多孔巖石中波致流體流動引起的頻散和衰減效應(yīng).這些數(shù)字巖心的尺度雖然比實際地層要小很多，但其揭示的基本物理過程和規(guī)律卻對綜合利用多頻段地球物理數(shù)據(jù)進行多尺度地質(zhì)特征刻畫有重要的啟示意義.

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

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

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

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

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

1 基本原理

1.1 邊界條件

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

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

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

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是質(zhì)點在x和z方向的速度分量，ux和uz是質(zhì)點在x和z方向的位移分量.

當(dāng)t≥0時，數(shù)字巖心左邊界質(zhì)點的速度分量為

vx(0,z,t)=0,

vz(0,z,t)=0.

(2)

當(dāng)t≥0時，數(shù)字巖心右邊界質(zhì)點的速度分量為

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

vz(L,z,t)=0.

(3)

式中，L表示數(shù)字巖心的邊長；f(t)是t的函數(shù)，f(0)=0；A是調(diào)節(jié)應(yīng)變大小的常數(shù)，使數(shù)字巖心在x方向滿足小應(yīng)變假設(shè).

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

圖2 歸一化的f(t)、f(t)的積分和f(t)的振幅譜(a) 函數(shù)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 波場模擬

假設(shè)數(shù)字巖心的骨架由各向同性的線彈性介質(zhì)組成，其本構(gòu)方程為

(4)

式中，λ和μ是拉梅常數(shù)；σxx，σzz和σzx表示應(yīng)力；exx，ezz和ezx表示應(yīng)變.

假設(shè)數(shù)字巖心飽和牛頓流體，牛頓流體的本構(gòu)關(guān)系為

(5)

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

用位移和應(yīng)力表示的運動方程為

(6)

式中，ρ表示密度.

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

圖3 旋轉(zhuǎn)交錯網(wǎng)格示意圖Fig.3 Sketch map of rotated staggered grid

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

(7)

式中，ψi可表示網(wǎng)格頂點的位移和速度，下標(biāo)i=1,2,3,4表示網(wǎng)格的4個頂點的標(biāo)號.

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

1.3 縱波的速度與逆品質(zhì)因子

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

(8)

1.4 精度驗證

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

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

圖4 牛頓流體的頻散和逆品質(zhì)因子(a) 速度頻散曲線； (b) 速度的數(shù)值解與解析解之差； (c) 逆品質(zhì)因子曲線； (d) 逆品質(zhì)因子的數(shù)值解與解析解之差.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 含裂隙致密巖石的頻散和衰減特征數(shù)值模擬

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

2.1 數(shù)字巖心構(gòu)建

為了刻畫擠噴流效應(yīng)，數(shù)字巖心包含若干相互獨立的裂隙—圓孔系統(tǒng).裂隙為直裂隙，走向垂直外部擠壓(拉伸)方向，一端與圓孔相連，另一端封閉，見圖5.數(shù)字巖心a和b的邊長為1.0×10-3m，孔隙度分別約為0.020和0.051.裂隙的長度為8.8×10-5m，寬度為3×10-6m，扁率(也稱縱橫比，定義為矩形裂隙的寬度與長度的比值)約為0.034，圓孔半徑為1.5×10-5m.數(shù)字巖心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).數(shù)字巖心a、b、c和d的裂隙密度分別約為0.035，0.089，0.048和0.104.

圖5 數(shù)字巖心數(shù)字巖心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.數(shù)字巖心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個數(shù)字巖心的骨架由同一固體介質(zhì)構(gòu)成，孔隙被同一流體飽和.固體與流體的參數(shù)見表1.綜合考慮模擬效果、裂隙扁率、模型大小和計算負(fù)擔(dān)，我們對流體設(shè)計了較大的黏滯系數(shù)(朱偉等，2020).

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

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

2.2 模擬測試

當(dāng)裂隙飽和黏性流體時，需要對裂隙寬度進行充分采樣以便以適當(dāng)?shù)木瓤坍嫴ㄖ铝餍?yīng).我們以數(shù)字巖心b為例進行網(wǎng)格剖分和數(shù)值模擬試驗，確定離散裂隙寬度所需的網(wǎng)格數(shù)量.首先，將裂隙寬度離散為3個網(wǎng)格；然后，以此為基礎(chǔ)，逐步加密網(wǎng)格(即減小網(wǎng)格大小)，分別得到裂隙寬度為6個和12個網(wǎng)格的模型.如此，我們得到3個離散的數(shù)字巖心網(wǎng)格模型，它們的大小和結(jié)構(gòu)完全相同，只是網(wǎng)格大小和數(shù)量不同.裂隙寬度分別為3個、6個和12個網(wǎng)格.

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

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

圖6 數(shù)字巖心b的頻散和逆品質(zhì)因子曲線曲線3-grids、6-grids和12-grids分別對應(yīng)裂隙寬度用3個、6個和12個網(wǎng)格離散.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.

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

與網(wǎng)格大小有關(guān)的另一不利影響因素是數(shù)值頻散.但當(dāng)波長與網(wǎng)格邊長的比值較大時，數(shù)值頻散可以勿略.當(dāng)數(shù)字巖心的裂隙寬度分別為3個，6個和12個網(wǎng)格時，網(wǎng)格邊長分別為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.波長與網(wǎng)格邊長的最小比值約為1517，即最短波長遠(yuǎn)遠(yuǎn)大于最大網(wǎng)格邊長，在此條件下數(shù)值頻散非常小，不影響擠噴流的頻散效應(yīng).

2.3 模擬結(jié)果

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

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

圖7 數(shù)字巖心的模擬結(jié)果(a) 速度; (b) 逆品質(zhì)因子.曲線a、b、c、d分別表示數(shù)字巖心a、b、c、d的模擬結(jié)果，φ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.

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

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

3 討論

3.1 單頻率與寬頻帶動態(tài)應(yīng)力應(yīng)變模擬的關(guān)系

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

圖8 單頻率和寬頻帶模擬結(jié)果比較(a) 速度; (b) 逆品質(zhì)因子.數(shù)字巖心為a和b，單頻率模擬的頻率為50 kHz、65 kHz和75 kHz.圖例中a-band和b-band表示數(shù)字巖心a和b的寬頻帶模擬結(jié)果，a-single和b-single表示數(shù)字a和b的單頻率模擬結(jié)果.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 含裂隙致密巖石的頻散和衰減特征的控制因素分析

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

(9)

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

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

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

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

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

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

4 結(jié)論

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