一維頻率域中觀尺度虛巖石物理方法

吳建魯 吳國忱

(①中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580;②海洋國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室，山東青島 266071)

1 引言

地震波在地下雙相介質(zhì)中傳播時會發(fā)生速度頻散和能量衰減。20世紀(jì)50年代以來，許多學(xué)者對地震波在雙相介質(zhì)中的衰減機理做了大量研究。從研究地震波衰減因素尺度的角度出發(fā)，可將地震波衰減機制分為三類，即宏觀尺度、微觀尺度和中觀尺度[1，2]。由宏觀尺度“Biot流”[3-5]和微觀尺度“擠噴流”[6]引起的地震波衰減都無法預(yù)測地震頻帶的衰減和頻散。White等[7，8]在考慮中觀尺度(遠大于孔隙尺寸且小于地震波長的尺度定義為中觀尺度)流體斑塊飽和情況時，預(yù)測了地震波在地震頻帶的衰減和頻散。隨后人們針對中觀尺度的衰減和頻散機理做了大量的理論研究[9-16]，結(jié)果表明地下介質(zhì)的中觀尺度非均質(zhì)性對地震波在地震頻帶的衰減起主控作用。

在構(gòu)建特定理想化巖石物理模型的基礎(chǔ)上發(fā)展特定的地震波衰減理論是研究地震波衰減和頻散特征的重要手段和方法。但是由于實際介質(zhì)的巖性、流體斑塊結(jié)構(gòu)、巖石孔滲性差異等因素的復(fù)雜變化，理論模型往往無法與實驗測量的結(jié)果進行精確的匹配分析，從而無法驗證理論的正確性和適用條件等[17-19]。Masson等[20，21]將低頻形式的動態(tài)滲透率關(guān)系式轉(zhuǎn)換至?xí)r間域后，采用準(zhǔn)靜態(tài)時間域有限差分交錯網(wǎng)格方法，利用Biot彈性波方程首次針對中觀尺度非均質(zhì)模型進行虛巖石物理理論分析，并將該方法推廣至全頻帶。針對中觀尺度非均質(zhì)性，Quintal等[22]和Rubino等[23]采用有限元法進行地震波衰減和速度頻散的研究。由于地下油氣儲層與地震頻帶的衰減和頻散有著密切的聯(lián)系，所以針對由中觀尺度非均質(zhì)性引起的地震波衰減和頻散研究、了解地震波在雙相介質(zhì)中的傳播規(guī)律十分必要[24-26]?；诘卣鸩ㄋp理論，前人也利用實際地震資料進行了地震波衰減信息提取及儲層預(yù)測方面的研究[27-30]。

本文借鑒準(zhǔn)靜態(tài)模擬思路，用頻率域有限差分法開展中觀尺度非均質(zhì)虛巖石物理模型的研究。與時間域虛巖石物理方法相比，該方法既可以直接求取任一頻率下的地震波衰減和頻散速度，便于應(yīng)用于實際巖樣的地震波衰減預(yù)測分析，又消除了針對頻率域的動態(tài)滲透率公式進行近似帶來的誤差，也避免了討論在巖樣外表面施加的力源函數(shù)表達式及時間剖分穩(wěn)定性等問題。首先通過數(shù)值計算模擬不同數(shù)量特征單元疊置模型的地震波衰減和頻散，說明了特征單元表征模型衰減信息的唯一性；隨后分析了相同含氣飽和度條件下不同非均勻尺度及斑塊結(jié)構(gòu)模型的地震波衰減和頻散；最后，計算了三相流體周期性層狀介質(zhì)的地震波衰減和頻散結(jié)果及模型的位移與應(yīng)力的頻率域空間分布特征，可知不同流體之間的中觀尺度相對流動是誘導(dǎo)地震波在地震頻帶衰減的主控因素。

2 一維形式雙相介質(zhì)理論

本文主要采用前人推導(dǎo)的兩種雙相介質(zhì)的控制方程：第一種為無耗散力的線性固結(jié)耦合的準(zhǔn)靜態(tài)方程[31]；第二種為考慮流體黏滯性導(dǎo)致能量耗散雙相介質(zhì)方程[3-5]。前者未考慮體力耗散，可以精確地描述地震頻帶內(nèi)波的傳播特征，后者考慮了由宏觀尺度流固相對運動摩擦導(dǎo)致的“Biot”流，進而產(chǎn)生地震波高頻衰減，前者是后者的低頻形式。通過對這兩種控制方程的研究，可以全面檢驗本文所提出的一維頻率域虛巖石物理方法的準(zhǔn)確性。

在假設(shè)孔隙介質(zhì)的應(yīng)力散度為零的前提下，線性固結(jié)耦合的準(zhǔn)靜態(tài)Biot方程在頻率—空間域的一維形式[32]為

(1)

Biot在考慮宏觀尺度“Biot流”的基礎(chǔ)上推導(dǎo)了飽和流體孔隙介質(zhì)條件下的本構(gòu)方程，初步建立了地震波在雙相介質(zhì)中傳播的理論框架，成為后人研究地震波在雙相介質(zhì)中傳播的理論基礎(chǔ)[2，3]。其含流體孔隙介質(zhì)彈性波傳播的一維控制方程為

(2a)

(2b)

(3)

3 一維頻率域準(zhǔn)靜態(tài)虛巖石物理方法

當(dāng)外界施加一定的力時，巖石樣本會產(chǎn)生形變(圖1)，通過對上述兩種孔隙介質(zhì)地震波控制方程微分項進行常規(guī)空間離散

(4)

(5)

式中L為模型的長度。利用壓力P0與體應(yīng)變θ的比值可求取模型的復(fù)變模量E，即

(6)

從而獲取模型的地震波衰減因子Q-1(ω)和頻散速度VP(ω)

(7)

在上、下施加垂直方向外力P0時，為保證巖石物理模型的封閉條件，使上、下邊界的流體位移Uf為零，則模型的上、下邊界條件為

(8)

與傳統(tǒng)的頻率域有限差分?jǐn)?shù)值模擬方法相比，上述的頻率域虛巖石物理方法不需要吸收邊界條件，減小了模型網(wǎng)格的剖分?jǐn)?shù)量，提高了數(shù)值計算效率。

圖1 一維虛巖石物理模型示意圖

4 數(shù)值計算結(jié)果

為驗證頻率域一維虛巖石物理方法的有效性，利用周期性層狀斑塊飽和模型(圖2a，介質(zhì)參數(shù)見表1)進行數(shù)值模擬，其非均質(zhì)尺度為中觀尺度。飽含兩種不同液相水平層狀孔隙介質(zhì)周期性交互排列，橫向無限延伸，僅考慮垂直方向一維應(yīng)力—應(yīng)變響應(yīng)。首先考慮含氣砂巖與含水砂巖周期排列情況，各層厚度均為0.1m，即含氣飽和度為0.5，空間網(wǎng)格步長dz設(shè)為1mm，采用該模型特征單元(圖2b)進行地震波衰減和頻散模擬，結(jié)果如圖3所示，其中曲線“Biot1941”和曲線“Biot1956”分別為基于準(zhǔn)靜態(tài)控制方程(式(1))和考慮宏觀尺度“Biot流”的雙相介質(zhì)控制方程(式(2))求取的結(jié)果，曲線“White1975”為利用White等[8]基于周期性層狀介質(zhì)的復(fù)變模量解析解的結(jié)果(解析解公式參見文獻[8，14])。由圖可見，地震波衰減和頻散對頻率有很強的依賴性，這與前人的實驗室測量結(jié)果[17]一致。與解析解相比，基于兩種雙相介質(zhì)傳播方程、采用本文所描述的頻率域有限差分方法可以獲取精確的地震波衰減和頻散結(jié)果，曲線“Biot1956”在高頻帶地震波速度迅速降低，因為在此頻帶內(nèi)地震波波長與模型尺度相當(dāng)，地震波產(chǎn)生共振現(xiàn)象，為地震波能量傳播停滯所致[6，10]。

表1 孔隙介質(zhì)參數(shù)

圖4為針對周期性層狀介質(zhì)斑塊飽和模型，基于準(zhǔn)靜態(tài)控制方程和考慮宏觀尺度“Biot流”的雙相介質(zhì)控制方程模擬不同頻率下特征單元(圖2b)的應(yīng)力、固體位移、流體位移及流體壓力的分布。由圖可見，基于準(zhǔn)靜態(tài)假設(shè)的雙相介質(zhì)控制方程正演模擬的應(yīng)力場為常數(shù)，與其方程描述的應(yīng)力場散度為零的條件一致，而基于考慮宏觀尺度“Biot流”的雙相介質(zhì)控制得到的應(yīng)力場在高頻振蕩變化，這直觀展示了地震波發(fā)生共振時的應(yīng)力場分布特征(圖4a)。當(dāng)某一特定頻帶內(nèi)固體位移、流體位移及流體壓力空間變化較大(即梯度較大)時，則產(chǎn)生的地震波衰減和頻散相對較強。本文僅討論地震頻帶附近的地震波衰減和頻散，即不考慮高頻帶宏觀尺度“Biot流”引起的地震波共振影響，故下文采用準(zhǔn)靜態(tài)Biot的雙相介質(zhì)控制方程。

分別對圖2b～圖2d描述的特征單元采用頻率域模擬方法求取地震波衰減和頻散，分析特征單元數(shù)量對周期性層狀介質(zhì)斑塊飽和模型的地震波衰減和頻散的影響(圖5)。由圖可見，周期性層狀介質(zhì)的斑塊飽和模型的地震波衰減和頻散不受特征單

圖2 周期性層狀斑塊飽和模型及其特征單元示意圖 (a)周期性層狀模型示意圖； (b)單個特征單元； (c)兩層特征單元疊置； (d)三層特征單元疊置

圖3 周期性層狀斑塊飽和介質(zhì)模型的地震波衰減(a)和頻散(b)模擬結(jié)果

圖4 基于準(zhǔn)靜態(tài)控制方程(左)和考慮宏觀尺度“Biot流”的雙相介質(zhì)控制方程(右)數(shù)值模擬結(jié)果 (a)應(yīng)力； (b)流體壓力； (c)固體位移； (d)流體位移

圖5 不同個數(shù)特征單元疊加的地震波衰減(a)和頻散(b)模擬結(jié)果

元數(shù)量的影響，驗證了本文數(shù)值求取方法的穩(wěn)定性，表明該模型表征的地震波衰減和頻散的唯一性。

圖6為相同含氣飽和度下不同斑塊分布的特征單元示意圖，含氣或含水的特征單元的總厚度設(shè)為0.12m，空間網(wǎng)格步長為1mm。特征單元一與周期性層狀介質(zhì)斑塊飽和模型的特征單元相同； 特征單元二上、下分別為厚度0.12m的含氣、含水孔隙介質(zhì)； 特征單元三、四則把含氣層與含水層分別進行二、四等分后交互排列，單個液相的介質(zhì)非均勻尺寸逐漸減??； 特征單元五、六、七、八中的含氣和水孔隙介質(zhì)分別按1∶2比例進行等分，然后進行如圖所示的排列組合。由于所構(gòu)建的上述特征單元不再是周期性層狀介質(zhì)模型所模擬的理想介質(zhì)，其解析解求取十分困難，所以采用本文提出的數(shù)值模擬算法求解各個特征單元的地震波衰減和頻散，結(jié)果如圖7所示。由特征單元一到特征單元四的地震波衰減和頻散結(jié)果可知： 當(dāng)非飽和流體的非均勻尺度逐漸減小時，地震波衰減主頻向高頻段移動，頻帶變寬；當(dāng)非飽和流體的非均勻尺度不變，流體斑塊尺度組合關(guān)系不同時，則求取的數(shù)值模擬結(jié)果會存在差異，如圖7中特征單元五～八的計算結(jié)果所示； 當(dāng)含氣飽和度及非均勻尺度不變時，其中下部含水孔隙介質(zhì)的層厚對地震波的衰減和頻散結(jié)果起主導(dǎo)作用，且由于存在多種形式非均質(zhì)組合，特征單元五、六的地震波衰減出現(xiàn)“多峰”現(xiàn)象。從不同特征單元的速度頻散結(jié)果(圖7b)可見，在相同含氣飽和度條件下，不同非均勻尺度情況地震波頻散速度曲線也存在差異。值得注意的是，不同特征單元頻散速度的高、低頻極限還是一致的。

圖6 相同含氣飽和度下不同斑塊飽和分布示意圖 淺色表示含氣層；深色表示含水層

圖7 不同特征單元的地震波衰減(a)和頻散(b)模擬曲線

為進一步討論多相孔隙介質(zhì)中觀尺度非均勻飽和情況的地震波衰減和頻散，采用頻率域模擬方法模擬含三相流體周期性層狀介質(zhì)模型(圖8)的地震波衰減和頻散及高、低頻極限(圖9)，并分析應(yīng)力場、固體位移、流體位移及流體壓力的變化情況[6](圖10)。模擬時油的參數(shù)設(shè)置為：Koil=1GPa；ρoil=700kg/m3；ηoil=0.04Pa·s。當(dāng)?shù)卣鸩ㄔ谠撃Ｐ椭袀鞑r，多相流體之間會存在不同情況中觀尺度“局域流”現(xiàn)象(圖10)，從而發(fā)生不同頻帶的地震

波衰減和頻散，即存在地震波衰減“雙峰”現(xiàn)象。高、低頻帶的地震波衰減分別是由氣相與油相、水相之間的相對流體產(chǎn)生的，由于油、水之間的物性(模量、黏滯性及密度等)差別不大，引起的地震波衰減較小，無法體現(xiàn)在計算結(jié)果之中，說明當(dāng)非飽和地下儲層賦存一定的天然氣時，儲層的地震波衰減特征比較明顯，有利于油氣儲層等方面的研究。

圖8 三相流體周期性層狀介質(zhì)示意圖

圖9 三相流體周期性層狀介質(zhì)的地震波衰減(a)和頻散(b)模擬曲線

圖10 三相流體周期性層狀介質(zhì)模型應(yīng)力場(a)、流體壓力(b)、固體位移(c)及流體相對位移(d)模擬結(jié)果

5 認識與討論

地震波在地下含流體孔隙介質(zhì)中傳播時，中觀尺度的“局域流”衰減機制是引起地震波在地震頻帶產(chǎn)生衰減和頻散的主要因素。在前人針對中觀尺度“局域流”的研究的基礎(chǔ)上，利用Biot(1941)準(zhǔn)靜態(tài)方程和考慮宏觀尺度“Biot流”的Biot(1956)雙相介質(zhì)方程，采用頻率域有限差分模擬方法，應(yīng)用周期性層狀斑塊飽和模型驗證了本文方法的有效性。本文方法不受模型特征單元數(shù)量影響，可以模擬任意形式斑塊分布的一維雙相介質(zhì)模型，突破了理論解析特征單元理想的假設(shè)條件。利用該方法進一步分析了相同含氣飽和度條件下不同結(jié)構(gòu)斑塊飽和尺度模型和三相流斑塊模型的地震波衰減和頻散。斑塊飽和非均質(zhì)結(jié)構(gòu)不同，產(chǎn)生的地震波衰減和頻散也不相同，當(dāng)非飽和流體儲層中含有一定的天然氣時，地震波衰減相對較大，表明了利用地震波衰減進行油氣預(yù)測的潛在價值。另一方面，該方法可以獲取巖石物理模型任意頻率的地震波衰減和頻散信息，便于與巖石物理實驗中測量結(jié)果進行對比分析。

本文針對一維模型的頻率域地震波衰減和頻散模擬方法精度較高，但模型為二維甚至三維時，由于地震波傳播方程的參數(shù)變量多，建立的阻抗矩陣大，存在一定的不適用性。這就需要減少波動方程的參數(shù)變量的個數(shù)，增加邊界條件的微分階數(shù)，以提高算法的適用性。

[1] Dvorkin J，Nur A.Dynamic poroelasticity：a unified model with the squirt and the Biot mechanisms.Geophysics，1993，58(4)：524-533.

[2] 張廣智，何鋒，張佳佳等.微觀與介觀波致流下的速度頻散與衰減.石油地球物理勘探， 2017，52(4)：743-751. Zhang Guangzhi，He Feng，Zhang Jiajia et al.Velocity dispersion and attenuation at microscopic and mesoscopic wave-induced fluid flow.OGP，2017，52(4)：743-751.

[3] Biot M A.Theory of propagation of elastic waves in a fluid-saturated porous solid Ⅰ：Low-frequency range.The Journal of the Acoustical Society of America，1956，28(2)：168-178.

[4] Biot M A.Theory of propagation of elastic waves in a fluid-saturated porous solid Ⅱ：Higher-frequency range.The Journal of the Acoustical Society of America，1956，28(2)：179-191.

[5] Biot M A.Mechanics of deformation and acoustic propagation in porous media.Journal of Applied Physics，1962，33(4)：1482-1498.

[6] Müller T，Gurevich B，Lebedev M.Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks — A review.Geophysics，2010，75(5)：A147-A164.

[7] White J E.Computed seismic speeds and attenuation in rocks with partial gas saturation.Geophysics，1975，40(2)：224-232.

[8] White J E，Mikhayloya N G，Lyakhovitskiy F M.Low-frequency seismic waves in fluid-saturated layered rocks.Physics of the Solid Earth，1975，11(10)：654-659.

[9] 席道瑛，邱文亮，程經(jīng)毅等.飽和多孔巖石的衰減與孔隙率和飽和度的關(guān)系.石油地球物理勘探， 1997， 32(2)：196-201. Xi Daoying，Qiu Wenliang，Cheng Jingyi et al.Relation between attenuation and porosity or saturation of rock.OGP，1997，32(2)：196-201.

[10] 劉炯，馬堅偉，楊慧珠.White球狀Patchy模型中縱波傳播研究.地球物理學(xué)報，2010，53(4)：954-962. Liu Jiong，Ma Jianwei，Yang Huizhu.Research on P-wave’s propagation in White’s sphere model with patchy saturation.Chinese Journal of Geophysics，2010，53(4)：954-962.

[11] 鄧?yán)^新，王尚旭，杜偉.介觀尺度孔隙流體流動作用對縱波傳播特征的影響研究——以周期性層狀孔隙介質(zhì)為例.地球物理學(xué)報，2012，55(8)：2716-2727. Deng Jixin，Wang Shangxu，Du Wei.A study of the influence of mesoscopic pore fluid flow on the propagation properties of compressional wave — a case of periodic layered porous media.Chinese Journal of Geo-physics，2012，55(8)：2716-2727.

[12] 吳國忱，吳建魯，宗兆云.周期性層狀含孔隙、裂隙介質(zhì)模型縱波衰減特征.地球物理學(xué)報，2014，54(8)：2666-2677. Wu Guochen，Wu Jianlu，Zong Zhaoyun.The attenuation of P wave in a periodic layered porous media containing cracks.Chinese Journal of Geophysics，2014，54(8)：2666-2677.

[13] Nicola T，Quintal B.Seismic attenuation in partially saturated rocks：Recent advances and future directions.The Leading Edge，2014，33(2)：640-646.

[14] 吳建魯，吳國忱，宗兆云.含混合裂隙、孔隙介質(zhì)的縱波衰減規(guī)律研究.地球物理學(xué)報，2015，57(4)：1378-1389. Wu Jianlu，Wu Guochen，Zong Zhaoyun.Attenuation of P waves in a porous medium containing various cracks.Chinese Journal of Geophysics，2015，57(4)：1378-1389.

[15] 李世凱，文曉濤，阮韻淇等.基于黏滯—彌散理論的含氣砂巖數(shù)值模擬與分析.石油地球物理勘探， 2017，52(4)：752-759. Li Shikai，Wen Xiaotao，Ruan Yunqi et al.Sandstone gas reservoir simulation and analysis based on diffusive-viscous theory.OGP，2017，52(4)：752-759.

[16] 陳程，文曉濤，郝亞炬等.基于White模型的砂巖儲層滲透率特性分析.石油地球物理勘探，2015，50(4)：723-729. Chen Cheng，Wen Xiaotao，Hao Yaju et al.Sandstone reservoir permeability characteristics analysis based on White model.OGP，2015，50(4)：723-729.

[17] Batzle M，Han D，Hofmann R.Fluid mobility and frequency-dependent seismic velocity — Direct measurements.Geophysics，2006，71(1)：N1-N9.

[18] Madonna C，Tisato N.A new seismic wave attenua-tion module to experimentally measure low-frequency attenuation in extensional mode.Geophysical Prospecting，2013，61(2)：302-314.

[19] McCann C，Sothcott J，Best A I.A new laboratory technique for determining the compressional wave properties of marine sediments at sonic frequencies and in situ pressures.Geophysical Prospecting，2014，62(1)：97-116.

[20] Masson Y J，Pride S R，Nihei K T.Finite difference modeling of Biot’s poroelastic equations at seismic frequencies.Journal of Geophysical Research：Solid Earth，2006，111(10)：B10305.

[21] Masson Y J，Pride S R.Poroelastic finite difference modeling of seismic attenuation and dispersion due to mesoscopic-scale heterogeneity.Journal of Geophysical Research：Solid Earth，2007，112(3)：B03204.

[22] Quintal B，Steeb H，F(xiàn)rehner M et al.Quasi-static finite element modeling of seismic attenuation and dispersion due to wave-induced fluid flow in poroelastic media.Journal of Geophysical Research：Solid Earth，2011，116(B1)：B01201.

[23] Rubino J G，Holliger K.Seismic attenuation and velocity dispersion in heterogeneous partially saturated porous rocks.Geophysical Journal International，2012，188(3)：1088-1102.

[24] 黃中玉，王于靜，蘇永昌.一種新的地震波衰減分析方法——預(yù)測油氣異常的有效工具.石油地球物理勘探，2000，35(6)：768-773. Huang Zhongyu，Wang Yujing，Su Yongchang.A new analysis method for seismic wave attenuation effective tool for predicting hydrocarbon anomaly.OGP，2000，35(6)：768-773.

[25] 張會星，何兵壽，姜效典等.利用地震波在雙相介質(zhì)中的衰減特性檢測油氣.石油地球物理勘探，2010，45(3)：343-349. Zhang Huixing，He Bingshou，Jiang Xiaodian et al.Utilizing attenuation characteristic of seismic wave in dual-phase medium to detect oil and gas.OGP，2010，

45(3)：343-349.

[26] 張大偉，孫贊東，王學(xué)軍等.利用零井源距VSP資料進行品質(zhì)因子反演.石油地球物理勘探，2011，46(S1)：47-52. Zhang Dawei，Sun Zandong，Wang Xuejun et al.Q-factor inversion and calculation using zero-offset VSP data.OGP，2011，46(S1)：47-52.

[27] 蘭鋒，李明，李艷東等.應(yīng)用黏彈性波動方程模擬地震低頻陰影.石油地球物理勘探，2015，50(1)：84-90. Lan Feng，Li Ming，Li Yandong et al.Low-frequency shadow numerical simulation with visco-elastic wave equation.OGP，2015，50(1)：84-90.

[28] Wu X，Chapman M，Li X et al.Quantitative gas saturation estimation by frequency-dependent amplitude-versus-offset analysis.Geophysical Prospecting，2014，62(6)：1224-1237.

[29] Innanen K A，Mahmoudian F.Characterizing the de-gree of amplitude-variation-with-offset nonlinearity in seismic physical modelling reflection data.Geophysical Prospecting，2014，62(6)：1-8.

[30] 胡軍輝，文曉濤，許艷秋等.利用黏滯—彌散波動方程理論進行油水識別.石油地球物理勘探，2016，51(3)：556-564. Hu Junhui，Wen Xiaotao，Xu Yanqiu et al.Oil-water recognition based on diffusive-viscous wave equation.OGP，2016，51(3)：556-564.

[31] Cavallini F，Carcione J M，Vidal De Ventós D et al.Low-frequency dispersion and attenuation in anisotropic partially saturated rocks.Geophysical Journal International，2017，209(3)：1572-1584.

[32] Biot M A.General theory of three-dimensional consoli-dation.Journal of Applied Physics，1941，12(2)：155-164.