基于數(shù)字巖心模擬的深層頁巖氣儲(chǔ)層可壓性評(píng)價(jià)模型

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

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

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

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ì)頁巖，有機(jī)質(zhì)豐富，巖心上還有筆石富集與少量黃鐵礦充填物。巖心整體頁理較發(fā)育，少量局部發(fā)育。微裂縫較發(fā)育且少量微裂縫呈高角度分布。發(fā)育有大量的微米—納米級(jí)無機(jī)孔隙和有機(jī)質(zhì)孔隙。儲(chǔ)層以石英和黏土礦物為主，兩者之和的平均體積分?jǐn)?shù)為40%，其中白云石、斜長(zhǎng)石、黃鐵礦、方解石等礦物僅占少量，地層非均質(zhì)性較強(qiáng)，石英含量由下往上逐漸減少。

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

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

多組分三維數(shù)字巖心是在構(gòu)建有機(jī)質(zhì)、黏土、方解石、石英、白云石、斜長(zhǎng)石、黃鐵礦單組分?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)建的多個(gè)不同4種成分的多組分?jǐn)?shù)字巖心，即可研究不同礦物組分對(duì)于巖石彈性參數(shù)的影響；另外，為了考慮層理與裂縫的影響，本文中建立平板狀的層理(圖2(k))、平板狀裂縫(圖2(m))巖心模型，將其嵌入多組分?jǐn)?shù)字巖心模型，從而建立層理發(fā)育(圖2(l))和裂縫發(fā)育(圖2(n)的三維數(shù)字巖心模型；為了考慮含油氣飽和度的影響，對(duì)于多尺度多組分巖心模型(圖2(j))中的孔隙空間(圖2(i))進(jìn)行氣水兩相流體驅(qū)替模擬獲取不同氣水分布巖心模型(圖2(p))，研究氣水飽和度對(duì)于巖石彈性參數(shù)的影響。

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

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

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

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

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

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

2.2 有限元法計(jì)算巖石彈性參數(shù)驗(yàn)證

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

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ù)字巖心有限元法模擬出來的彈性模量值可靠，且理論計(jì)算彈性模量與有限元計(jì)算彈性模量的平均相對(duì)誤差為：體積模量5.55%，剪切模量2.60%。

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

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

圖6 試驗(yàn)測(cè)試值與有限元計(jì)算值的平均相對(duì)誤差Fig.6 Average relative error between experimental test value and finite element calculation value

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

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

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

(4)

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

(5)

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

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

(6)

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

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

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

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

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 可壓性評(píng)價(jià)模型建立

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

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

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

根據(jù)前文的影響因素敏感性分析，確定本文中所研究的9種影響因素對(duì)于彈性參數(shù)變化敏感性的大小，因此根據(jù)其對(duì)可壓性的相對(duì)重要性進(jìn)行取值可得到可壓性指標(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，再對(duì)特征向量作歸一化處理，便可得到各參數(shù)的權(quán)重系數(shù)為0.312、0.222、0.156、0.108、0.074、0.051、0.035、0.025、0.018。故最終可建立可壓性評(píng)價(jià)模型。裂縫傾角、層理傾角、裂縫密度、裂縫長(zhǎng)度是通過對(duì)成像測(cè)井中的裂縫和層理進(jìn)行人機(jī)交互識(shí)別得到的，裂縫密度為單位長(zhǎng)度內(nèi)裂縫發(fā)育的條數(shù)，層理密度通過成像測(cè)井紋層識(shí)別后用紋層密度來表征，紋層密度為單位長(zhǎng)度內(nèi)紋層的個(gè)數(shù)；孔隙半徑從核磁共振測(cè)井中獲得，主要利用其橫向弛豫時(shí)間(T2譜)來反映巖石的孔徑分布，即巖石在飽和流體狀態(tài)下T2譜越長(zhǎng)則巖石孔隙度越大；脆性礦物含量通過巖心試驗(yàn)分析得到?？蓧盒栽u(píng)價(jià)模型為

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為歸一化裂縫長(zhǎng)度；Kw為歸一化脆性礦物；Kh為歸一化水平應(yīng)力差異系數(shù)；Cl,dip為歸一化層理傾角；Kj為歸一化孔隙半徑；Bhd為歸一化含氣飽和度。

4 可壓性測(cè)井評(píng)價(jià)應(yīng)用

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

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

由于成像測(cè)井的分辨率很大，所以對(duì)可壓性特征明顯的部分做了截取展示(圖9(b)、圖10(b))。從計(jì)算的可壓性曲線FI可知，X1井龍-14小層到五峰組整體的可壓性指數(shù)都較高，范圍為0.42～0.69，平均為0.58，龍-11小層到五峰組頂部最高(圖9)；而X3井在龍-11小層和五峰組的可壓性指數(shù)相對(duì)較高(X3井沒有核磁測(cè)井?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儲(chǔ)層可壓性特征相對(duì)較好。

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

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

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

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

5 結(jié) 論

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

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