基于孔隙縱橫比譜反演的飽和巖石寬頻段巖石物理模型

韓旭，王尚旭，劉浩杰，劉韜，唐跟陽*

1 中國石油大學(北京) 地球物理學院，北京 102249

2 中國石化勝利油田物探研究院，山東 257000

3 中國石化石油勘探開發研究院，北京 100728

0 引言

巖石的速度主要受孔隙度、礦物成分、流體類型和飽和度控制[1]。大量實驗室研究表明，干巖石的彈性速度或彈性模量隨著有效壓力的增加而顯著增加[2-4]。這種壓力依賴性通常可以通過裂縫根據其縱橫比逐漸閉合來解釋[4-5]。雖然這些裂縫只占孔隙空間的一小部分，但相比于球形孔，它們對多孔巖石的有效速度有更加顯著的影響。通常來說，Gassmann理論模型由于忽略巖石孔隙結構影響，往往不能解釋巖石彈性性質的壓力依賴性。在實驗室條件下使用超聲波測量的實驗數據證實了孔隙結構對聲速的重要性[6-8]。已經建立了許多理論公式來將有效彈性特性與孔隙結構相關聯。由于巖石中孔隙空間的復雜性，通常假設它們可以由以其縱橫比為特征的橢球體表示，這些橢球狀孔隙可以表征各種各樣的孔隙形狀，并且可以通過Eshelby(1957)的單一橢球體等效夾雜理論進行分析處理[9]。典型的此類理論包括 K-T 模型和 DEM 理論[4,10-12]，它們將橢圓形孔隙視為巖石固相中的包含物。Shapiro(2003)利用雙重孔隙結構模型解釋巖石速度隨壓力的這種非線性/線性變化現象，利用孔隙縱橫比α表征孔隙形狀[13]。孔隙可以視為具有不同形狀的硬孔隙(α＞ 0.01)和軟孔隙(α＜0.01)的組合[14]。這種雙重孔隙度模型很好地解釋了有效壓力下速度的非線性到線性變化的特征。Cheng and Toks?z (1979)提出了一種與壓力相關聯的孔隙縱橫比和孔隙閉合的關系式，通過超聲測量的干燥和飽和速度反演孔隙和孔隙縱橫比的離散型分布[15]。該模型適用于超聲頻段，并且不考慮流體效應。

另一方面，對于流體飽和巖石，壓力和頻率的影響是耦合的，因為它們通過孔隙的微觀結構相互連接[16]。不同形狀和方向的孔隙之間的局部流體流動是流體飽和巖石中彈性波的頻散和衰減的主要原因之一[17]。孔隙結構對巖石彈性參數的影響會因流體與裂孔微觀結構的相互作用而進一步復雜化。在噴射流模型中，軟孔孔隙度是造成彈性速度頻散的主要因素之一，噴射流特征頻率與軟孔隙的縱橫比有關[18]。Gurevich et al(2010)提出了一種簡單噴射流模型[19]，其中孔隙的縱橫比以及微裂縫的含量是造成飽和流體巖石頻散與衰減的主要因素。但是，該理論模型只是給出了單一孔隙縱橫比以及微裂隙含量對彈性模量頻散與衰減的影響，該假設不符合真實巖石的孔隙結構。并且，這些理論都有許多表征孔隙形狀的“自由”參數(例如孔隙縱橫比、軟孔孔隙度等)。這些參數是未知的，通常需要擬合實驗數據獲得的。這會降低這些模型的擬合能力和實用性，因為良好的擬合不足以證明模型的適用性。

本文的目標是克服這種不確定性，并從速度或彈性模量的壓力依賴性獲得此類參數。通常，孔隙形狀和孔隙度可以表征巖石彈性參數的壓力依賴性，可以使用對這種依賴性的分析來估計這些參數。Shapiro (2003)提出的干速度壓力依賴性理論模型[13]。然而，模型中的軟孔隙假設具有相同的縱橫比。這種將巖石的孔隙縱橫比簡化為單一縱橫比的方法無法完全表征巖石微裂隙的分布，通常巖石的孔隙縱橫比具有較寬的范圍。利用實驗測量獲得的巖石速度與壓力的相關性，并計算巖石軟孔隙縱橫比的分布，可以更為有效的獲取巖石孔隙縱橫比的分布。Cheng and Toks?z(1979)將實驗室測量的各種巖石的地震速度作為壓力和飽和度相關的函數進行反演，得到它們的孔隙形狀譜[15]。這種反演是基于孔隙閉合作為壓力和速度變化的函數的理論。通過超聲測量的干燥和飽和速度反演孔隙度和孔隙縱橫比的離散型分布，但是該方法忽略了流體流動效應。

本文擴展了原有的巖石物理模型。首先，利用超聲測量的速度-壓力變化曲線來反演孔隙縱橫比分布及其孔隙度。其次，基于K-T等效介質模型將干燥硬孔隙加入巖石基質中，并再次加入飽和軟孔隙(軟孔隙中添加流體頻變體積模量)獲得巖石“干骨架”的彈性模量；最后利用Gassmann流體替換理論描述硬孔隙中飽和流體對彈性模量的影響。

1 巖石物理實驗

為了研究復雜孔隙結構巖石速度的壓力和頻率依賴性，本文選取了中國東部油田沙河街組的致密砂巖樣品進行超聲速度測量。致密砂巖樣品的孔隙度為4%，密度為2644 kg/m3。樣品的礦物含量通過XRD衍射儀獲得，并利用VRH等效介質理論模型計算樣品的基質體積模量(Ks=55 GPa)和剪切模量(μs=29.6 GPa)。樣品所飽和的液體為68#白油(體積模量2 GPa，密度890 kg/m3，黏度0.068 Pa·s)。采用超聲透射法測量不同壓力下不同流體飽和樣品的彈性參數。

超聲波透射法是利用一對超聲波換能器激發和接收在巖石樣品中傳播的超聲波P波或S波，記錄接收到的波形并用于速度計算。 換能器的直徑為17 mm，激勵頻率約為1 MHz。在測量過程中，溫度保持在20 ℃，壓力從0 MPa增加到50 MPa。首先，對致密砂巖樣品進行潤濕以軟化“過干燥”樣品。將其放入充滿水蒸氣的密封容器中48小時，以軟化巖石骨架。在下一步中，使用“偽自吸”方法將樣品浸入油中，該方法首先將樣品浸入油中，抽真空，然后在50 MPa下加壓長達7天。然后，測量了油飽和致密砂巖的超聲波速度。

圖1顯示了不同壓力下致密砂巖樣品的縱波和橫波速度。比較了這些測量結果，發現在干燥條件下，超聲速度的變化隨壓力從非線性增加變為幾乎線性增加，轉變壓力約為30 MPa。當樣品充滿液體時，速度隨著壓力的增加而變化更加平滑，不會表現出非線性/線性轉變。隨著壓力的增加，干燥和油飽和條件下的縱波速度之間的差異減小，在干燥和油飽和條件下橫波速度差小于縱波速度差。速度曲線表明壓力對超聲速度的影響很大，隨著壓力的增加，巖石樣品中的軟孔隙逐漸閉合[14]。在低壓范圍(0～30 MPa)的非線性趨勢下，縱橫比小的軟孔隙逐漸閉合，表明軟孔隙度是速度變化的主要控制因素。在高壓范圍(30～50 MPa)的線性趨勢下，軟孔隙幾乎閉合，表明硬孔隙是速度變化的主要控制因素[14,16]。

2 巖石物理建模

2.1 孔隙縱橫比譜反演

上述致密砂巖實驗結果表明，速度隨壓力的非線性到線性變化趨勢與雙重孔隙結構(球型硬孔隙與裂縫型軟孔隙)模型一致。為此，遵循 Cheng and Toks?z (1979) 提出的反演方法[15]，從超聲波速度-壓力曲線關系中反演孔隙縱橫比譜。對于具有巖石基質(巖石礦物顆粒)模量Ks，μs和孔隙填充物模量Ki，μi的兩相材料，Kuster and Toks?z (1974)給出了等效模量K*和μ*的表達式：

其中，αm是第m個縱橫比，c(αm)是它在相關的孔隙度。T1和T2是由背景材料的體積模量Ks和剪切模量μs以及孔隙填充物的體積模量Ki和剪切模量μi組成的表達式，該參數反映孔隙形狀影響的耦合因子。在 Cheng and Toks?z (1979)中給出了它們的明確表達[15]。

根據Eshelby (1957)對橢圓形包含物的定義[9]，Toks?z 等(1976)給出了孔隙體積分數dc/c隨壓力P變化的表達式[20]。

其中，Ei是縱橫比和有效模量的表達式，KA是靜態干體積模量，通常用超聲波速度計算的干體積模量近似。在Cheng and Toks?z (1979)中給出了它們的明確表達[15]。通常認為d/ 1c c＜－，假設裂縫是閉合的。微裂隙在壓力Pn下的大小與其在零壓力下的大小的關系為：

對于縱橫比為α和半長軸為a的扁球狀孔隙，c(α) =4παa3/3。孔隙縱橫比隨壓力的變化由下式給出：

因此，可以將孔隙縱橫比的變化與其在零壓力下的狀態聯系起來：

將方程(4)代入方程(1)和(2)并得到：

等式(7)和(8)具有y=Ax的形式，其中，Kn*和μn*是超聲條件在不同壓力下測量的彈性參數。為了使矩陣A更加穩定，選擇參數向量x為c(αm)/αm而不是c(αm)。這樣，模型向量x的元素大約在0.1和1之間。方程總數為2LN，其中L表示飽和度數，N表示壓力數。m是壓力點的數量或縱橫比。模型向量x的維數為m-1，因為它應滿足所有縱橫比的孔隙度之和等于總孔隙度ctotal。

系數矩陣A的維數為2LN×m-1，因此系統y=Ax是一個超定線性方程組。通過最小二乘法轉換A[21]。

其中，δ是與ATA的最小特征值相關聯的可調參數。

圖2a顯示了致密砂巖樣品在0 MPa下的離散孔隙縱橫比譜。縱橫比的分布顯示了典型的雙孔隙模型的特征。硬孔隙在致密砂巖的總孔隙度中占主導地位。在零壓力下，軟孔隙度具有相當寬的縱橫比范圍(10-3～10-4.5)，軟孔孔隙度為0.03%，而硬孔孔隙度為3.98%。圖2b為不同壓力下的孔隙縱橫比譜。硬孔孔隙度隨壓力的增加幾乎保持不變，而軟孔孔隙度的數量隨著壓力的增加逐漸減少，數值逐漸減小。表明這些軟孔隙隨著壓力的增加而逐漸閉合(孔隙度小于10-8基本閉合)。圖3為不同壓力條件下，軟、硬孔隙度的絕對值變化曲線。如圖3所示，硬孔隙基本不隨壓力的變化而變化，而軟孔隙度隨著壓力的增加逐漸降低，當壓力達到50 MPa，軟孔隙度基本為0。

圖2 致密砂巖樣品反演孔隙縱橫比譜結果Fig.2 Inversion pore aspect ratio spectrum results of tight sandstone samples

圖3 軟、硬孔孔隙度隨壓力的變化曲線Fig.3 Variation of compliant and stiff porosity with pressure

圖4是干燥條件下超聲實驗測試結果以及K-T模型的模擬結果。基于離散孔隙縱橫比譜的K-T模型對速度有很好的擬合結果。由于線性迭代反演方法依賴于有效介質理論，因此沒有考慮流體流動效應。如果除了干數據之外還使用飽和巖石樣品的速度，它會在孔隙縱橫比譜的反演中引入誤差。利用干燥數據反演的孔隙縱橫比譜能夠更加準確的計算巖石的孔隙結構。當樣品只含有硬孔隙時，樣品的縱橫波速度不隨這壓力的變化而變化，在壓力范圍(40～50 MPa)模擬結果和實驗結果有很好的擬合關系。然而，在低壓力范圍，模擬結果與實驗結果差異較大，因此，軟孔隙度是巖石速度的主要控制因素。

圖4 基于干燥數據反演孔隙縱橫比譜的K-T模型模擬結果，紅線為只含有硬孔隙的模擬結果Fig.4 K-T model simulation results of inversion of pore aspect ratio spectrum based on drying data, the red line is the simulation result with only stiff pores

2.2 “干骨架”彈性模量計算

對于飽和流體巖石，流體流動機制是解釋巖石速度頻散的主要理論之一。基于K-T等效介質模型和反演的孔隙縱橫比譜的結果，計算巖石樣品的“干骨架”彈性模量。其中，硬孔隙中不添加流體，軟孔隙中加入流體頻變體積模量[22]。

對于液體飽和巖石，超聲頻率下的修正骨架模量與干骨架模量相關[19]，公式為：

其中，Kmf和μmf是修正干骨架的體積模量和剪切模量。Kh是不含有軟孔隙巖石的體積模量，cc是軟孔孔隙度。

對于中間頻率，Gurevich 等(2010) 表明模量遵循相同的方程(11)-(12)但流體模量Ki由與頻率相關的(ω)代 替[19]。Glubokovskikh等(2016)簡 化 和 概 括了Gurevich et al.(2010) 的方法[22]。對于具有動態黏度η的孔隙流體Kf*(ω)，其相應表達式由Tsai and Lee (1998) 給出[23]：

其中，Jk是k階的第一類修正貝塞爾函數，而μf*(ω) =iωη。

首先，基于K-T模型將硬孔隙加入到巖石基質中，硬孔隙的含量可通過反演的孔隙縱橫比譜獲得。

其中，Kh和μh是只含有硬孔隙巖石的彈性模量，α13－是模擬出的硬孔隙縱橫比，c(α13－)是硬孔孔隙度。這里我們只把縱橫比大(α＞0.01)的硬孔孔隙度帶入到K-T模型中，可獲得不含有微裂隙巖石的干骨架模量Kh和μh。

其次，在高頻極限下，修正的“干骨架”的體積和剪切模量由等式(1)-(2)定義。繼續添加飽和流體頻變體積模量的軟孔隙，通過用Kf*(ω)和μf*(ω)替換方程(1)-(2)中的流體體積模量Ki和μi來獲得較低頻率的模量，其中，Kh和μh可以通過方程(15)-(16)計算得到。

2.3 全飽和巖石彈性模量計算

最后，全飽和巖石彈性模量由Gassmann流體替代方程計算[24]。

式中，Ksat和μsat是全飽和巖石體積模量和剪切模量。cs是硬孔孔隙度。

圖5a是新模型和K-T模型對超聲速度的擬合結果。需要說明的是，孔隙縱橫比譜的反演不涉及全飽和速度數據。在超聲頻段，K-T模型和噴射流模型對P波和S波都具有很好的擬合結果。圖5(b)是噴射流模型、K-T模型和Gassmann模型對體積模量頻散的模擬結果。如圖所示，在超聲頻率范圍內，流體飽和樣品的建模速度對于新模型和等效介質模型(K-T)沒有太大差異。在低頻范圍內，流體飽和樣品的Gassmann模型與改進模型的低頻極限相同。因此，新模型的高頻極限符合等效介質模型，而低頻極限符合Gassmann理論。相比于簡單噴射流模型，新模型由于考慮了孔隙縱橫比譜，頻散范圍更寬。

圖5 (a)致密砂巖超聲速度與改進K-T模型和K-T模型的擬合結果；(b) 4種模型模擬致密砂巖體積模量隨頻率的變化關系Fig.5 (a) The ultrasonic velocity of the tight sandstone and the fitting results of the improved K-T model and the K-T model; (b) The relationship between the bulk modulus of the tight sandstone simulated by the four models with frequency

圖6顯示了不同壓力條件下，P波速度和S波速度頻散的結果。106Hz的數據點是對于壓力和油飽和條件下的超聲P波和S波速度。通過孔隙縱橫比分布，新模型可用于模擬彈性參數的壓力和頻率依賴性。從模擬結果可以看出，速度的特征頻率隨壓力增加而減小，而速度頻散量降低。隨著壓力的增加，樣品的軟孔隙的平均縱橫比和孔隙度不斷下降。較低的縱橫比會導致特征頻率降低，而軟孔孔隙度的降低會導致頻散幅度下降[19,25-28]。

圖6 不同壓力條件下速度頻散的模擬結果 (a)縱波速度; (b)橫波速度Fig.6 Simulation results of velocity dispersion under different pressure conditions (a) P-wave velocity; (b) S-wave velocity

3 討論

基于不同壓力下的孔隙縱橫比譜，在干燥和完全飽和條件下，新模型可以很好地預測致密砂巖P波和S波速度(圖4a)。在超聲頻率范圍內，對于新建立的巖石物理模型和有效介質模型，流體飽和樣品的建模速度沒有太大差異。這里同樣采用Gassmann理論模型來解釋致密砂巖飽和油的超聲速度。結果表明， Gassmann理論模型忽略了孔隙結構的影響，模擬結果完全低估了實際測量結果。因此，相比于Gassmann理論模型，新建立的模型具有更高的精度。同時，該模型的一個優點是不需要可調參數。輸入參數和變量來自獨立的實驗室測量，建模結果具有可驗證性。實際測量和建模結果都表明，壓力和頻率對彈性模量的影響是耦合的，因為這兩個因素與孔隙微觀結構是緊密聯系在一起的。壓力的變化將引起孔隙結構的改變(孔隙縱橫比和軟孔孔隙度)，而孔隙結構的改變會影響巖石速度頻散的大小和特征頻段。另外由于缺乏低頻數據，可以利用低頻數據進一步驗證模型的適用性。

4 結論

巖石微觀孔隙結構不僅是干燥巖石速度的主要控制因素，也決定了飽和流體巖石中與流體流動相關的速度頻散和衰減作用。為了研究壓力和頻率對流體飽和巖石彈性模量的影響，本文根據干燥條件下超聲速度-壓力曲線變化來反演孔隙縱橫比譜。基于這種復雜的孔隙結構，本文擴展了原有的模型來解釋速度的壓力和頻率依賴性。在所提出的模型中，孔隙空間由兩部分組成：具有縱橫比分布的軟孔隙和縱橫比大的硬孔隙。通過定量表征軟孔隙閉合過程，獲得了軟孔隙每個壓力下的縱橫比分布。在軟孔隙中加入了流體頻變體積模量，使模型具有頻散效應。相比于簡單噴射流模型，新模型由于考慮了孔隙縱橫比譜，頻散范圍更寬。模型對與壓力相關的超聲數據有很好的擬合結果。新建立的模型對速度頻散給出合理的理論解釋。因此，該研究有助于我們更好地了解流體誘發地震波衰減和頻散的特征。