賴富強，劉粵蛟，張海杰，王海濤，毛海艷，張國統(tǒng)

(1.重慶科技學(xué)院復(fù)雜油田勘探開發(fā)重慶市重點實驗室，重慶 401331；2.重慶頁巖氣勘探開發(fā)有限責(zé)任公司，重慶 401147)

深層頁巖氣儲層一般是指深度大于3 500 m的頁巖氣儲層[1]，由于沉積作用、成巖作用、埋深、溫度及壓力等因素影響，深層頁巖氣儲層相較中淺層層理、斷裂及天然微裂縫更為發(fā)育，應(yīng)力差及巖石力學(xué)參數(shù)變化更大，孔隙結(jié)構(gòu)更加復(fù)雜，導(dǎo)致儲層的可壓性主控因素不清[2-4]。而且由于頁巖儲層為多尺度、多組構(gòu)儲集空間，具有非均質(zhì)性強、易風(fēng)化等特點[5-7]，受限于巖心來源、收獲率等因素，難以取得具有代表性的樣品，開展傳統(tǒng)巖石物理實驗的難度較大，所以需要構(gòu)建深層頁巖氣儲層典型的多尺度、多組分三維數(shù)字巖心模型來進行巖石參數(shù)模擬[8-10]。國內(nèi)外學(xué)者從巖石力學(xué)性質(zhì)、礦物組分、孔隙結(jié)構(gòu)、層理及裂縫發(fā)育和地質(zhì)因素等方面對儲層進行可壓性評價分析[11-12]，但是由于儲層各異，每種模型考慮的可壓性影響因素存在差別，而且并未給出如何確定每種影響因素對于可壓性變化的敏感性大小，主控因素權(quán)重不清楚，只是定性地規(guī)定了各影響因素之間的重要程度，所以構(gòu)建的可壓裂性評價模型的準(zhǔn)確性和適用性較差[12]。筆者提出綜合考慮研究區(qū)深層頁巖氣儲層特點，由于沒有頁巖巖心三維X射線CT掃描，利用已有的二維巖心圖片，利用馬爾科夫鏈-蒙特卡洛法[14](MCMC法)構(gòu)建多尺度、多組分三維數(shù)字巖心，基于數(shù)字巖心物理模型，對深層頁巖氣儲層可壓性主控因素采用有限元計算方法進行彈性參數(shù)數(shù)值模擬[15]，然后利用摩爾斯分類篩選法定量地分析裂縫傾角、裂縫密度、層理密度、裂縫長度、脆性礦物、層理傾角、孔隙半徑、含氣飽和度這8種影響因素引起彈性參數(shù)變化的敏感性相對大小。最后綜合考慮水平應(yīng)力差異系數(shù)的影響，利用層次分析法確定這9種影響因素(上述8種影響因素和水平應(yīng)力差異系數(shù))各自的權(quán)重，建立一種深層頁巖儲層的可壓性評價計算模型。

1 區(qū)域概況

研究區(qū)為大足區(qū)塊，位于渝西地區(qū)北部，地處四川盆地丘陵地帶，東南部為低山，最高海拔934 m，中部為淺丘帶壩，西部為深丘，最低海拔為185 m[16]。四川盆地大地構(gòu)造分區(qū)包括川西北坳陷帶、川中隆起帶和川東南坳褶帶(圖1)，大足區(qū)塊在構(gòu)造上屬于川東南坳褶帶、川中隆起帶交界[17]。區(qū)內(nèi)穩(wěn)定發(fā)育有埋深介于3 500～4 500 m的五峰組-龍馬溪組頁巖層，是頁巖氣開發(fā)的主要目標(biāo)層序。

圖1 區(qū)域構(gòu)造位置示意圖Fig.1 Schematic diagram of regional structure location

根據(jù)研究區(qū)域X1井、X2井、X3井、X5井的資料分析，龍馬溪組發(fā)育有深灰色、黑灰色、灰黑色及黑色頁巖，自上而下顏色逐漸變深，主要為黑色碳質(zhì)頁巖，有機質(zhì)豐富，巖心上還有筆石富集與少量黃鐵礦充填物。巖心整體頁理較發(fā)育，少量局部發(fā)育。微裂縫較發(fā)育且少量微裂縫呈高角度分布。發(fā)育有大量的微米—納米級無機孔隙和有機質(zhì)孔隙。儲層以石英和黏土礦物為主，兩者之和的平均體積分?jǐn)?shù)為40%，其中白云石、斜長石、黃鐵礦、方解石等礦物僅占少量，地層非均質(zhì)性較強，石英含量由下往上逐漸減少。

2 頁巖氣儲層可壓性模擬與主控因素敏感性

通過對研究區(qū)大量的測井資料對比分析認(rèn)為，研究區(qū)可壓性評價的關(guān)鍵參數(shù)包括了礦物組分、水平地應(yīng)力差、層理特征、裂縫、孔隙結(jié)構(gòu)，含油氣飽和度等[11-13]。在沒有巖心三維CT掃描信息的情況下，本文中針對研究區(qū)巖心鑄體薄片(圖2(a))進行二維掃描電鏡與X射線衍射實驗，對薄片礦物組分(石英、黏土、方解石、白云石、斜長石、黃鐵礦)及有機質(zhì)進行識別得到二維圖片(圖2(b))，對二維圖像進行遍歷掃描獲得每個方向上的條件概率，然后依據(jù)MCMC法對3個相互垂直的二維圖片(圖2(c))進行三維重構(gòu)[18]單組分?jǐn)?shù)字巖心(圖2(d)～(g))。

多組分三維數(shù)字巖心是在構(gòu)建有機質(zhì)、黏土、方解石、石英、白云石、斜長石、黃鐵礦單組分?jǐn)?shù)字巖心的基礎(chǔ)上通過選取其中4種單組分?jǐn)?shù)字巖心嵌套得來(圖2(h))。多尺度多組分三維數(shù)字巖心(圖2(j))是將高分辨率孔隙信息(圖2(i))融合到低分辨率多組分三維數(shù)字巖心得來[19]。本文中所構(gòu)建的數(shù)字巖心的分辨率為3 μm/像素，模型尺寸為240 μm×240 μm×240 μm。

圖2 重構(gòu)的多尺度多組分三維數(shù)字巖心Fig.2 Reconstructed multi-scale and multi-component 3D digital core

基于構(gòu)建的多個不同4種成分的多組分?jǐn)?shù)字巖心，即可研究不同礦物組分對于巖石彈性參數(shù)的影響；另外，為了考慮層理與裂縫的影響，本文中建立平板狀的層理(圖2(k))、平板狀裂縫(圖2(m))巖心模型，將其嵌入多組分?jǐn)?shù)字巖心模型，從而建立層理發(fā)育(圖2(l))和裂縫發(fā)育(圖2(n)的三維數(shù)字巖心模型；為了考慮含油氣飽和度的影響，對于多尺度多組分巖心模型(圖2(j))中的孔隙空間(圖2(i))進行氣水兩相流體驅(qū)替模擬獲取不同氣水分布巖心模型(圖2(p))，研究氣水飽和度對于巖石彈性參數(shù)的影響。

用MCMC法構(gòu)建數(shù)字巖心的優(yōu)勢在于其根據(jù)原始圖像的轉(zhuǎn)移概率對新圖像進行賦值重構(gòu)，間接地考慮了巖心內(nèi)部物質(zhì)的空間分布，可以較好地表征實際巖心。

通過計算本文中所構(gòu)建的數(shù)字巖心與二維圖像的自相關(guān)函數(shù)，分析數(shù)字巖心模擬的平均孔徑與巖心CT實驗平均孔徑的相對誤差，驗證了本文MCMC法構(gòu)建的數(shù)字巖心可以較好地反映研究區(qū)實際頁巖的孔隙結(jié)構(gòu)和礦物組成(圖3)。在建立可靠的數(shù)字巖心模型的基礎(chǔ)上，利用有限元計算方法計算巖心彈性參數(shù)(圖2(q))，采用摩爾斯分類篩選法計算各種可壓性影響因素變化引起的彈性參數(shù)變化的相對大小(圖2(r))，從而計算敏感性因子，確定各影響因素敏感性相對大小，為層次分析法確定權(quán)重及可壓性評價模型的構(gòu)建做準(zhǔn)備(圖2(s))。

圖3 自相關(guān)分析與平均孔徑對比Fig.3 Comparison between autocorrelation analysis and average aperture size

2.1 有限元法計算巖石彈性參數(shù)

給定數(shù)字巖心模擬時所考慮物質(zhì)的彈性參數(shù)——體積模量和剪切模量[20-23]，依據(jù)有限元法即可計算出巖石體積模量和剪切模量[24]。有限元計算彈性參數(shù)的原理是：巖心離散化，對每個網(wǎng)格進行彈性參數(shù)賦值；在巖心X、Y、Z方向上分別施加一個宏觀應(yīng)力，依據(jù)應(yīng)力應(yīng)變關(guān)系確定每個單元格的機械能，然后依據(jù)變分原理即能量取極小值，依據(jù)共軛梯度法求解獲取每個單元格上的位移分布，從而可以確定巖心總的位移，最終確定巖心應(yīng)變。已知應(yīng)力和應(yīng)變，便可確定巖心的彈性體積模量和剪切模量(圖4)。

圖4 有限元法計算巖石彈性模量原理Fig.4 Principle of calculating rock elastic modulus by finite element method

2.2 有限元法計算巖石彈性參數(shù)驗證

首先利用理論值進行驗證，將有限元計算的彈性模量，與V-R-H(Voight-Reuss-Hill)理論計算的等效彈性模量[25](理論值)進行對比。

Voight上限值：

(1)

Reuss下限值：

(2)

V-R-H值：

MVRH=(MV+MR)/2 .

(3)

式中,N為巖石中的組分?jǐn)?shù)；fi為第i種組分的體積分?jǐn)?shù)；Mi為第i種組分的體積模量或剪切模量，GPa。

從圖5中可以看出，模擬值均落在理論值上下限值之間，說明通過數(shù)字巖心有限元法模擬出來的彈性模量值可靠，且理論計算彈性模量與有限元計算彈性模量的平均相對誤差為：體積模量5.55%，剪切模量2.60%。

圖5 理論計算彈性模量與有限元計算彈性模量Fig.5 Elastic modulus by theoretical calculation and finite element calculation

然后用試驗參數(shù)值進行驗證，用研究區(qū)儲層4 037～4 038 m深度的12塊巖心的多組分巖心模擬結(jié)果與對應(yīng)深度巖心的三軸巖石力學(xué)參數(shù)試驗測的彈性模量、泊松比以及計算的脆性指數(shù)進行對比(圖6)。結(jié)果顯示三軸巖石力學(xué)參數(shù)實驗測的彈性參數(shù)與有限元計算結(jié)果的平均相對誤差為:彈性模量8.9%，泊松比5.6%，脆性指數(shù)7.8%。通過理論值與試驗值的驗證，說明基于構(gòu)建三維數(shù)字巖心物理模型，利用有限元法計算巖石彈性參數(shù)的結(jié)果可靠，可開展下一步彈性參數(shù)影響因素敏感性分析工作。

圖6 試驗測試值與有限元計算值的平均相對誤差Fig.6 Average relative error between experimental test value and finite element calculation value

2.3 頁巖彈性參數(shù)影響因素敏感性

由數(shù)字巖心巖石物理模擬出體積模量K與剪切模量G。由下式計算出彈性模量E和泊松比μ:

E=9KG/(3K+G),

(4)

μ=(3K-2G)/(6K+2G).

(5)

采用摩爾斯分類篩選法[26]，定量模擬影響因素在不同數(shù)值情況下彈性參數(shù)的變化大小定量分析影響因素的敏感性，彈性模量的模擬如圖7所示。

采用自變量以固定步長變化，靈敏度判別因子取多個平均值，公式如下：

(6)

式中,S為靈敏度判別因子；n為模型運行次數(shù)；Yi、Yi+1分別為模型第i次、第i+1次運行輸出值；Y0為計算結(jié)果初始值；pi、pi+1分別為第i次、第i+1次模型運算參數(shù)值相對初始參數(shù)值的變化百分率。

通過層理、脆性礦物、裂縫以及孔隙含氣飽和度對于巖心彈性參數(shù)的模擬及敏感性分析結(jié)果(圖7、8)可知：8種影響因素對于彈性參數(shù)的敏感性由大到小依次是裂縫傾角、裂縫密度、層理密度、裂縫長度、脆性礦物、層理傾角、孔隙半徑、含氣飽和度。

圖7 影響因素變化引起彈性模量變化的模擬Fig.7 Simulation of Youngs modulus change caused by change of influencing factors

另外，在頁巖儲層實際壓裂施工中裂縫在縱向上受最小水平主應(yīng)力的影響而變化并沿著最大水平主應(yīng)力方向進行延伸，即水平主應(yīng)力差較小時，人工裂縫會更多的與天然裂縫溝通，表現(xiàn)出較好的可壓裂性，因此本文中還把應(yīng)力因素同上面8種因素進行考慮，但是由于應(yīng)力很難用數(shù)字巖心來模擬，本文中對其敏感性的相對大小進行了9次排列，然后計算水平應(yīng)力差異系數(shù)在不同排列情況下的可壓性指數(shù)。經(jīng)微地震監(jiān)測和壓裂施工曲線驗證，水平應(yīng)力差異系數(shù)排在脆性礦物與層理傾角之間時與微地震監(jiān)測和壓裂施工曲線的對應(yīng)效果最好，故本文中定性地將水平應(yīng)力差異系數(shù)對可壓性的敏感性定在脆性礦物與層理傾角之間。水平應(yīng)力差異系數(shù)計算公式為

Kh=(σH-σh)/σh.

(7)

式中，Kh為水平差應(yīng)力系數(shù)；σH為水平最大主應(yīng)力，MPa；σh為水平最小主應(yīng)力，MPa。

圖8 X3井巖心彈性參數(shù)敏感性Fig.8 Sensitivity of core elastic parameters of well X3

3 可壓性評價模型建立

采用層次分析法綜合裂縫傾角、裂縫密度、層理密度、裂縫長度、脆性礦物、水平地應(yīng)力差、層理傾角、孔隙半徑、油氣飽和度建立一種能夠?qū)搸r氣儲層可壓性進行全面科學(xué)評價的計算模型。為得到綜合評價儲層可壓性的無量綱常數(shù)，還要對各參數(shù)采用經(jīng)驗賦值或極差變換來進行標(biāo)準(zhǔn)化、歸一化處理。

(1)可壓性評價參數(shù)標(biāo)準(zhǔn)化?？蓧盒灾笖?shù)的影響指標(biāo)包括正向指標(biāo)、逆向指標(biāo)，正向指標(biāo)的值越大，逆向指標(biāo)的值越小，可壓裂性越好。正向指標(biāo)Si和逆向指標(biāo)Sj分別為

Si=(Xi-Xmin)/(Xmax-Xmin),

(8)

Sj=(Xmax-Xj)/(Xmax-Xmin).

(9)

式中，Xmax為參數(shù)最大值;Xmin為參數(shù)最小值；Xi、Xj為考慮參數(shù)的參數(shù)值。

(2)因素權(quán)重的確定及計算。利用層次分析法，根據(jù)參數(shù)敏感性的大小，確定它們之間重要性(表1)，并根據(jù)其重要性給每一層元素給出定量表征，構(gòu)造出判斷矩陣A，Aij表示元素i相對于元素j的重要程度值，Aji表示元素j相對于元素i的重要程度值。然后運用數(shù)學(xué)的方法求出判斷矩陣的最大特征值所對應(yīng)的特征向量,從而得到權(quán)重系數(shù)。

表1 判斷矩陣元素的標(biāo)度極其含義Table 1 Scale and meaning of judgment matrix elements

根據(jù)前文的影響因素敏感性分析，確定本文中所研究的9種影響因素對于彈性參數(shù)變化敏感性的大小，因此根據(jù)其對可壓性的相對重要性進行取值可得到可壓性指標(biāo)的判斷矩陣，如表2所示。

表2 判斷矩陣元素的標(biāo)度極其含義Table 2 Scale and meaning of judgment matrix elements

求得判斷矩陣的特征向量為W=(0.710 9，0.506 4，0.354 1，0.244 9，0.168 2，0.115 4，0.079 7，0.056 3，0.041 8)，最大特征值λmax=9.401 4，再對特征向量作歸一化處理，便可得到各參數(shù)的權(quán)重系數(shù)為0.312、0.222、0.156、0.108、0.074、0.051、0.035、0.025、0.018。故最終可建立可壓性評價模型。裂縫傾角、層理傾角、裂縫密度、裂縫長度是通過對成像測井中的裂縫和層理進行人機交互識別得到的，裂縫密度為單位長度內(nèi)裂縫發(fā)育的條數(shù)，層理密度通過成像測井紋層識別后用紋層密度來表征，紋層密度為單位長度內(nèi)紋層的個數(shù)；孔隙半徑從核磁共振測井中獲得，主要利用其橫向弛豫時間(T2譜)來反映巖石的孔徑分布，即巖石在飽和流體狀態(tài)下T2譜越長則巖石孔隙度越大；脆性礦物含量通過巖心試驗分析得到。可壓性評價模型為

FI=0.312Lf,dip+0.222Lfd+0.156Cld+0.108LfL+

0.074Kw+0.051Kh+0.035Cl,dip+0.025Kj+0.018Bhd.

(10)

式中，Lf,dip為歸一化裂縫傾角；Lfd為歸一化裂縫密度；Cld為歸一化層理密度；LfL為歸一化裂縫長度；Kw為歸一化脆性礦物；Kh為歸一化水平應(yīng)力差異系數(shù)；Cl,dip為歸一化層理傾角；Kj為歸一化孔隙半徑；Bhd為歸一化含氣飽和度。

4 可壓性測井評價應(yīng)用

根據(jù)以上對可壓性的研究內(nèi)容，利用本文中所建立的可壓性計算模型對X1、X3井進行可壓性系數(shù)計算，結(jié)果如圖9、10所示。

圖9 X1井可壓性指數(shù)計算結(jié)果Fig.9 Calculation diagram of compressibility index of well X1

由于成像測井的分辨率很大，所以對可壓性特征明顯的部分做了截取展示(圖9(b)、圖10(b))。從計算的可壓性曲線FI可知，X1井龍-14小層到五峰組整體的可壓性指數(shù)都較高，范圍為0.42～0.69，平均為0.58，龍-11小層到五峰組頂部最高(圖9)；而X3井在龍-11小層和五峰組的可壓性指數(shù)相對較高(X3井沒有核磁測井?dāng)?shù)據(jù)，沒有考慮孔徑參數(shù)的影響)，范圍為0.42～0.69，平均為0.58(圖10)。兩口井在龍-11小層和五峰組的平均可壓性指數(shù)為0.56。根據(jù)Rickman等[27]的可壓性判斷標(biāo)準(zhǔn)，綜合考慮X1、X3井的可壓性指數(shù)，認(rèn)為研究區(qū)五峰組-龍-11儲層可壓性特征相對較好。

圖10 X3井可壓性指數(shù)計算結(jié)果Fig.10 Calculation diagram of compressibility index of well X3

利用本文模型計算X3井直改平井段的可壓性指數(shù)，首先發(fā)現(xiàn)7、11、14段的可壓性指數(shù)較高，5、21段的可壓性指數(shù)較低。從可壓性指數(shù)較高的第7段的加沙壓裂施工曲線來看(圖11(a)、(c))，第7段180 min前有多處裂縫破裂顯示(有破裂顯示：排量不變，油壓上升又下降；油壓迅速下降，排量上升；泵壓不變，排量上升。180 min后砂堵進行沖砂，不做分析)；而可壓性指數(shù)較低的第21段并沒有什么裂縫破裂顯示(無破裂顯示：泵壓隨排量的上升而上升)[28]，油壓曲線和排量曲線都沒有變化(圖11(b)、(c))。各段所對應(yīng)的微地震監(jiān)測SRV(頁巖儲層改造體積)值也表現(xiàn)出同樣的高低特征(圖12)。說明本文模型所計算的可壓性指數(shù)可以與實際生產(chǎn)所測參數(shù)對應(yīng)，模型具有一定的可靠性，可以較好地為頁巖氣儲層射孔提供參考并促進頁巖氣開發(fā)。

圖11 X3井第7段及第21段壓裂施工曲線Fig.11 Fracturing operation curve of section 7th and 21th of well X3

圖12 X3井(直改平)儲層綜合分析Fig.12 Comprehensive analysis of reservoir in well X3 (straight to flat)

5 結(jié) 論

(1)通過數(shù)字巖心模擬，可以定量地分析彈性參數(shù)各影響因素的敏感性相對大小，可以為層次分析法建立可壓性評價模型時更準(zhǔn)確地確定各參數(shù)之間的權(quán)重。

(2)經(jīng)微地震監(jiān)測和壓裂施工曲線驗證，基于數(shù)字巖心模擬的可壓性評價模型可為深層頁巖氣儲層靶點層位選取和壓裂工藝優(yōu)化提供指導(dǎo)。但是本次研究在數(shù)字巖心模擬層理與裂縫的時候是通過構(gòu)建一個同分辨率的平板來實現(xiàn)，較為簡單，不能準(zhǔn)確地反應(yīng)真實的巖心層理、裂縫發(fā)育情況，需要后期研究進一步改進。