富有機質巖石橫波速度預測方法

劉致水 劉俊州 董 寧 包乾宗 王震宇 時 磊

(①長安大學地球物理系，陜西西安 710054； ②中國石化石油勘探開發研究院，北京 100083)

0 引言

富有機質頁巖、泥巖是近年來油氣勘探開發的熱點目標，常規地震資料解釋、AVO反演等仍是富有機質巖石“甜點”預測中應用最廣泛的技術。作為井震結合、疊前反演、脆性因子計算及應力評價的重要基礎數據之一，精確的橫波速度對有效識別頁巖“甜點”具有重要作用[1]。目前，針對砂巖、碳酸鹽巖儲層的巖石物理橫波速度預測方法[2-5]較為成熟，而針對富有機質巖石的相關研究較少。富有機質巖石因含有一定的有機質(干酪根)成分而區別于砂巖、碳酸鹽巖儲層。此外，部分富有機質巖石表現為低孔隙度，如Bakken頁巖[6]；也有一些表現為中—高孔隙度，如Monterey頁巖[6]。有機質的彈性特征與流體近似，與礦物差別較大，例如:鹽水的體積模量為2.2GPa、剪切模量為0、密度為1.03g/cm3，黏土的體積模量是21GPa、剪切模量是9GPa、密度是2.55g/cm3，而干酪根的體積模量為2.9GPa、剪切模量為2.7GPa、密度為1.30g/cm3[7]。雖然干酪根與流體的體積模量、剪切模量和密度數值差別較小，但是兩者有本質不同，即干酪根是具有剪切性的固體，而流體不具有剪切性[1]。因此，在構建富有機質巖石的物理模型過程中，需謹慎考慮干酪根的描述方式。

巖石中的孔隙形態十分復雜[8-9]，對巖石速度影響極大：球形孔隙(如印模孔)使巖石速度變高；扁狀孔隙(如裂縫)使巖石速度變低[2,4]。據此認識，學者針對砂巖、碳酸鹽巖構建了多個描述孔隙形態、孔隙度與速度之間關系的函數(巖石物理模型)，并以此為基礎進行橫波速度預測[10-13]。在針對富有機質巖石的研究過程中，學者們參考砂巖、碳酸鹽巖的巖石物理建模思路進行了諸多有益的探索[14-18]。對富有機質巖石的電子顯微鏡( BSE-SEM)掃描觀察顯示，干酪根以斑塊狀分布于巖石礦物顆粒間,其形態多種多樣[19-21]，這種特征使巖石的彈性參數分散。Sayers[14]認為每個富有機質巖石樣品中的干酪根形態不同，假設干酪根為硬幣形狀包含物，利用硬幣的縱橫比表征巖石中干酪根的形態變化，該方法能有效解釋具有低—超低孔隙度、高干酪根含量特征的Bakken頁巖的彈性參數分散現象。但由于只考慮干酪根及其形態的作用，忽略了孔隙及其形狀對巖石速度的影響，因此上述方法在低—超低孔隙度巖石中應用效果較好，在孔隙度稍高的巖石中則無法應用。董寧等[15]利用SCA-DEM模型把孔隙加入巖石基質中，再利用固體替代理論將干酪根加入含孔隙巖石中，這種方法考慮了孔隙形狀，但是沒有考慮干酪根顆粒的形態，且計算過程中將干酪根與孔隙分步加入巖石中，因此忽略了干酪根與孔隙之間的相互作用。Guo等[16]將干酪根和孔隙看作包含物，基于SCA巖石物理模型將干酪根與礦物、孔隙結合構建富有機質巖石的物理模型，該方法將干酪根等效為一種固定形態的顆粒，將孔隙假設為多種形態，雖然考慮了孔隙形態的復雜性，卻沒有考慮不同干酪根顆粒之間形態差異。與Guo等[16]的方法類似，Liu等[18]將干酪根作為固定形態的包含物，利用DEM模型將干酪根與礦物、孔隙結合構建富有機質巖石的物理模型，通過試算、統計誤差給定整個井段研究目標的固定干酪根顆粒形態參數。在富有機質巖石的物理模型構建及速度預測研究中，如何有效兼顧孔隙形狀及干酪根形態因素的影響，關于這個問題學界研究較少。

本文基于Kuster-Toks?z(KT)模型構建了一種富有機質巖石的橫波速度預測方法，該方法利用硬幣形狀包含物的縱橫比表征干酪根形態和孔隙形狀，在縱波速度約束下同時求算巖石中的等效干酪根顆粒和等效孔隙縱橫比，在反演參數基礎上預測橫波速度。因針對每個樣點求算得到干酪根顆粒形態和孔隙形狀，提高了巖石物理模型對富有機質巖石的描述精度，從而降低了橫波速度預測誤差。需要說明的是，本文方法是預測垂直于頁巖層理的速度，在直井中可與目前富有機質巖石“甜點”預測中廣泛應用的地震資料解釋、疊前疊后反演等各向同性方法匹配較好，在斜井中則需要先進行斜井的速度各向異性校正[22]。

1 基本原理

1.1 KT巖石物理模型

等效介質巖石物理模型，如KT模型[8]、微分等效介質理論[23]、自洽模型[7]等，認為巖石的彈性模量不僅與巖石基質礦物、包含物的種類和含量有關，還與包含物的幾何形態有關。其中，KT模型給出了實驗室高頻條件下，基質礦物中包含多種類型包含物時等效介質的體積模量和剪切模量的計算公式[7-8]

(1)

(2)

硬幣形狀包含物的形狀因子為[24]

(3)

(4)

在求算低頻條件下飽和流體巖石的彈性模量時，一般使用KT模型求算干巖石的彈性模量，再用Gassmann方程[25]計算飽和流體巖石的彈性模量。

1.2 富有機質巖石的等效介質模型

將富有機質巖石等效為由巖石基質、干酪根顆粒、孔隙(干孔隙或飽和流體孔隙)組成的混合物(圖1)。把干酪根和孔隙分別等效為隨機分布的、具有單一縱橫比的硬幣形狀固體顆粒和孔隙，其形態由硬幣形狀包含物的縱橫比表征：當縱橫比接近1時，顆粒形態接近球形；當縱橫比接近0時，顆粒形態為裂縫形。設等效干酪根顆粒縱橫比為αk、等效孔隙縱橫比為αp，N=2，則式(1)、式(2)可寫為

圖1 富有機質巖石的等效介質模型示意圖

=Vk(Kk-Km)Pmk+φ(Kf-Km)Pmf

(5)

=Vk(Gk-Gm)Qmk+φ(-Gm)Qmf

(6)

式中：Vk指干酪根的體積含量；φ指孔隙度；Kk、Gk指干酪根的體積模量和剪切模量；Kf指流體的體積模量；Pmk、Qmk指等效干酪根顆粒的形狀因子，是關于αk的函數；Pmf、Qmf指等效孔隙的形狀因子，是關于αp的函數。

得到彈性模量后，可求取巖石的縱、橫波速度

(7)

(8)

式中巖石密度ρ=ρkVk+ρfφ+ρm(1-Vk-φ)，其中ρk、ρf、ρm分別是干酪根、流體、礦物基質的密度。

圖2是利用式(5)～式(8)得出的低頻條件下飽和鹽水巖石的縱、橫波速度隨αk、αp變化的規律。該算例參考Monterey頁巖的孔隙度與干酪根含量。假設φ=0.08，Vk=0.08，背景基質為泥質，其體積模量、剪切模量分別取Km=39.54GPa、Gm=25.68GPa，密度取ρm=2.64g/cm3；干酪根顆粒的體積模量、剪切模量分別取Kk=2.9GPa、Gk=2.7GPa，密度取ρk=1.30g/cm3；假設孔隙中含鹽水，其體積模量Kf=2.65GPa，密度ρf=0.99g/cm3，令αk和αp均在0.001～1.000之間變化。圖2顯示等效干酪根顆粒縱橫比和等效孔隙縱橫比與巖石速度之間的關系是非線性的。抽取圖2中的部分數據并在圖3中顯示。由圖可見，在相同干酪根含量和孔隙度的情況下，等效干酪根顆粒和等效孔隙縱橫比都會對速度造成較大影響，且縱橫比越接近于1，速度越大；等效孔隙縱橫比對速度的影響程度大于等效干酪根縱橫比，原因在于干酪根的體積模量大于鹽水，以及干酪根的剪切模量不為0。

圖2 縱(a)、橫(b)波速度隨等效干酪根顆粒縱橫比αk、等效孔隙縱橫比αp變化曲面

圖3 縱(a)、橫(b)波速度隨等效干酪根顆粒縱橫比αk、等效孔隙縱橫比αp變化曲線抽取圖2中幾組數據

1.3 等效干酪根和孔隙縱橫比反演

式(5)和式(6)表明，富有機質巖石的彈性模量不僅是礦物基質、干酪根、孔隙流體體積分數及其彈性模量的函數，還與巖石的等效干酪根顆粒縱橫比αk以及等效孔隙縱橫比αp密切相關。將式(5)、式(6)代入式(7)、式(8)中可以建立vP、vS與αk、αp的非線性關系式[vP,vS]=f(αk，αp)，利用此式，既可根據巖石組分體積含量、彈性參數、αk、αp正演計算巖石的縱、橫波速度，也可以根據巖石的縱波速度(或縱、橫波速度聯合)反演αk、αp。利用式(5)～式(8)反演求取αk、αp時的目標函數為

OF=Wp|(vPt-vPc)|/vPt+

Ws|(vSt-vSc)|/vSt

(9)

該式是關于αk、αp的二元非線性函數。式中：vPt、vSt為測量的縱、橫波速度；vPc、vSc為預測的縱、橫波速度；WP、WS為加權因子，且滿足WP+WS=1.0。當只有縱波資料時，WP=1.0，WS=0；同時有縱、橫波資料時，可取WP=WS=0.5。

使用Vernik等[6]的一個實測數據點為例說明求解目標函數OF的過程。該樣點的實測孔隙度φ=0.043，Vk=0.182，縱、橫波速度分別為vPt=3.70km/s、vSt=2.43km/s，密度為ρ=2.43g/cm3；背景基質和干酪根的參數與圖2一致，令αk和αp都在0.001～1.000之間變化。

圖4為WP=1.0和WP=0.5時的兩種目標函數曲線，其中，vPc和vSc根據式(5)～式(8)求得。圖4顯示在定義域內目標函數表現為二維曲面，極值點位于藍色凹槽內。將圖4所示的數據在αk、αp兩個方向上取極小值并在圖5顯示，可見當WP=1.0和WP=0.5時，目標函數OF的形態基本一致；OF在0.001～1.000內的局部極值點[αk，αp]可取多個；OF的全局極值點大約為[αk，αp]=[0.035，0.070]，說明利用縱波速度約束求解目標函數與利用縱、橫波速度聯合約束求解目標函數的結果接近，也從側面說明利用縱波速度求算αk、αp，進而預測橫波速度的方法是可行的。本文采用一種非線性全局尋優粒子群算法——混沌量子粒子群算法[26]求解式(9)。當缺乏橫波測井資料時，可由縱波反演求取αk、αp，再將αk、αp代入式(5)～式(8)計算橫波速度；反之，利用縱、橫波聯合反演求得的αk、αp計算縱、橫波速度，并可以用來評估測量速度的質量。由于實際測量的縱、橫波速度數據常存在噪聲，特別是在井孔擴徑層段，聲波時差曲線出現周波跳躍，橫波速度受擴徑的影響更大，這種情況下實測數據不準確，需要進行必要的校正[27]。

圖4 反演目標函數隨αk、αp的變化曲面(a)縱波速度約束反演(Wp=1.0,Ws=0)目標函數； (b)縱、橫波速度聯合約束反演(Wp=0.5,Ws=0.5)目標函數

圖5 縱波速度約束與縱、橫波速度聯合約束目標函數隨αp(a)和αk(b)的變化曲線對比

2 實驗室數據測試

2.1 數據介紹

使用Vernik等[6]在高頻條件測量到的三組富有機質巖石數據(圖6)進行測試。這些數據在常溫、有效壓力70MPa條件下測得，樣品為干巖石。由圖6可見，Bakken頁巖的孔隙度為0.0104～0.0197，干酪根含量為0.122～0.423，屬低孔、高干酪根含量頁巖；Bazhenov頁巖的孔隙度和干酪根含量分別為0.0199～0.0420、0.066～0.207，屬低孔、高干酪根含量頁巖，但其孔隙度高于Bakken頁巖，干酪根含量低于Bakken頁巖；Monterey頁巖的孔隙度和干酪根含量分別為0.043～0.309、0.016～0.363，屬高孔頁巖，干酪根含量變化較大。

圖6 實測Bakken頁巖、Bazhenov頁巖、Monterey頁巖的孔隙度與干酪根含量交會圖

2.2 速度預測結果

利用本文方法對上述三組樣品進行試算。圖7為使用縱波速度約束和縱、橫波速度聯合約束下得到的預測與實測橫波速度交會圖。由圖可見，數據均勻分布于對角線(圖中紅色線)附近，說明預測結果與實測結果吻合度較高。

圖7 預測與實測橫波速度交會圖(a)縱波速度約束； (b)縱、橫波速度聯合約束

使用預測與實測速度之間的相對誤差平均值MAE、均方根誤差RMSE、相關系數R2等三項指標定量評價預測結果的可靠性

(10)

(11)

(12)

表1給出了利用縱波速度約束和利用縱、橫波速度聯合約束所得橫波速度的誤差統計結果，可以看到：①針對三組巖樣，利用縱波速度約束求算的結果與實測數據之間誤差小，說明利用縱波資料反演αk和αp，進而預測橫波速度的方法是可行的；②縱、橫波速度聯合約束反演所得結果的統計誤差優于縱波單獨約束，其原因在于橫波速度的參與。

表1 縱波速度約束與縱、橫波速度聯合約束所得橫波速度與實測橫波速度的誤差統計

圖8a為根據縱波速度約束和縱、橫波速度聯合約束反演的αk與干酪根含量的交會圖，圖8b為反演的αp與孔隙度的交會圖。結果表明：利用縱波速度約束反演和利用縱、橫波速度聯合約束反演所得結果較為接近，說明單獨利用縱波速度約束反演αp和αk是可行的。此外，對縱、橫波速度聯合約束反演的結果進行分析發現：①Bakken頁巖的αk均小于0.2，而αp均大于0.4，說明Bakken頁巖的等效干酪根顆粒以接近扁平的近裂縫形態為主，而其孔隙則接近球形；②Bazhenov巖樣中，4個樣品的αk在0.1以下，4個樣品的αk為0.4～0.7，6個樣品的αp小于0.1，2個樣品的αp大于0.9，說明其等效干酪根顆粒和孔隙同時存在接近裂縫的形態和接近球形的形態；③Monterey巖樣的αk為0～1.0，αp為0～0.8，說明其等效干酪根顆粒和孔隙也存在同時接近裂縫的形態和接近球形的形態；④對三個區塊的頁巖整體考慮，可以發現一個總趨勢：干酪根含量高的頁巖αk小；孔隙度大的頁巖αp小。

圖8 根據縱波速度約束和縱、橫波速度聯合約束反演的等效干酪根顆粒縱橫比與干酪根含量的交會圖(a)及等效孔隙縱橫比與孔隙度的交會圖(b)

2.3 對比方法及其計算結果

將本文方法與文獻中的三種單一參數自適應方法進行對比，以說明同時反演αk、αp在橫波速度預測中的作用。三種單一參數方法說明如下。

方法1將干酪根作為基質礦物的一種(不考慮干酪根顆粒的形態)，將孔隙等效為硬幣形狀且令孔隙縱橫比自適應變化。計算流程簡述為：①通過Voigt-Reuss-Hill(VRH) 平均公式[7]將基質礦物與干酪根混合；②利用KT模型將具有可變縱橫比的孔隙加入巖石基質中，在縱波速度約束下求算孔隙縱橫比；③將計算得到的孔隙縱橫比代入KT模型，計算富有機質巖石的彈性模量，進而求算速度。

方法2與本文流程近似。將干酪根與孔隙等效為硬幣狀顆粒，令孔隙縱橫比固定為常數(通過在有橫波資料的井中做實驗標定確定)，令干酪根縱橫比自適應。計算流程簡述為：①利用KT模型將具常數縱橫比的孔隙和具可變縱橫比的干酪根顆粒加入巖石基質中，在縱波速度約束下求算等效干酪根縱橫比；②將計算得到的等效干酪根縱橫比和常數孔隙縱橫比代入KT模型，計算富有機質巖石的彈性模量，進而求算速度。該方法賦予所有巖樣以相同的干酪根顆粒縱橫比，沒有考慮不同樣點間干酪根顆粒可能存在差別。

方法3與方法2流程近似，不同點在于令干酪根縱橫比固定為常數，而令孔隙縱橫比自適應。其計算流程為：①利用KT模型將具有常數縱橫比的等效干酪根顆粒和具有可變縱橫比的孔隙加入巖石基質中，在縱波速度約束下求算等效孔隙縱橫比；②將計算得到的等效孔隙縱橫比和常數等效干酪根縱橫比代入KT模型，計算等效介質巖石彈性模量，進而求算速度。

這種方法賦予所有巖樣以相同的孔隙縱橫比，沒有考慮不同樣品之間孔隙可能存在差別。

利用圖6所示數據進行試算。在此例中，對于方法2和方法3，使每組巖樣的固定參數在0～1之間以0.05的步長變化，在縱波速度的約束下求取最優可變參數。圖9為利用方法2、方法3預測的橫波速度誤差統計參數隨αk或αp的變化，由圖可見:對于方法2，當孔隙縱橫比分別為0.05、0.35、0.10時，通過干酪根縱橫比的自適應可使三組頁巖的預測橫波速度誤差分別取得最優值；對于方法3，當干酪根縱橫比分別為0.10、0.05、0.05時，通過孔隙縱橫比的自適應可使得三組頁巖的預測橫波速度的誤差分別取得最優值。說明要精準求取速度，不同頁巖選取的參數是不同的，這也證明了同時反演兩個參數進行橫波速度預測的必要性。

表2為三種單一參數自適應方法和本文方法預測橫波速度的誤差對比，其中方法2和方法3的αk或αp取圖9中的最優值。可以看出，相比三種單一參數自適應方法，本文方法的統計參數有明顯改善。

表2 四種方法預測橫波速度與實測速度的誤差統計

圖9 利用方法2(左)、方法3(右)預測的橫波速度與實測速度的誤差統計參數隨αk或αp的變化(a)MRE； (b)RMSE； (c)R2

3 實際測井數據測試

將本文方法應用到湖北省建南構造頁巖氣區塊A井侏羅系東岳廟段。圖10展示了目的層段的巖性錄井、儲層位置、測井解釋結果、速度預測結果及孔隙與干酪根縱橫比計算結果。利用多礦物測井解釋方法[28]估算孔隙度、礦物含量和含氣飽和度，解釋所得礦物主要為石英、長石、方解石、黏土。計算過程中使用的巖石組分的彈性參數和密度如表3所示，利用CARBOLOG方法[29]求算干酪根的體積分數，利用VRH平均公式[7]計算巖石基質的體積模量和剪切模量，利用Wood方程[7]計算混合流體的體積模量。將礦物、干酪根的體積含量、孔隙度、流體飽和度數據引入巖石物理模型，在縱波速度的約束下反演αk和αp，進而預測橫波速度。

表3 巖石組成成分的彈性參數和密度[1，7]

由圖10可見，四種方法預測的結果與實測數據趨勢基本一致，而新方法所得結果吻合度最高。表4為四種方法在實際資料中應用所得結果誤差統計，表明新方法預測結果優于其他三種對比方法。

圖10 建南構造頁巖氣區域A井侏羅系東岳廟段富有機質泥頁巖速度預測值與實測值對比

表4 四種方法應用于實際資料時的誤差統計

4 結束語

本文提出一種根據縱波速度同時反演富有機質巖石中等效干酪根顆粒縱橫比和等效孔隙縱橫比，進而開展橫波速度預測的方法。該方法利用KT模型建立礦物、干酪根、孔隙、流體、干酪根顆粒與孔隙縱橫比之間的數學關系，使用非線性全局尋優算法反演等效干酪根顆粒縱橫比和等效孔隙縱橫比，將反演得到的兩種縱橫比代入巖石物理模型中計算橫波速度。通過實驗室與實際測井資料的應用，證明本文方法預測的結果與實測數據吻合度較高。將該方法預測結果與其他三種單一自適應參數方法進行對比，結果顯示本文方法預測橫波速度的統計參數具有明顯優勢，證明了該方法的可行性及可靠性。