應(yīng)用核磁T2譜與數(shù)字巖心技術(shù)計算粒度分布方法

孫建孟,趙建鵬,閆偉超,劉學(xué)鋒

(1.中國石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,山東青島 266580；2.中國石油大學(xué)理學(xué)院,山東青島 266580)

應(yīng)用核磁T2譜與數(shù)字巖心技術(shù)計算粒度分布方法

孫建孟1,趙建鵬1,閆偉超1,劉學(xué)鋒2

(1.中國石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,山東青島 266580；2.中國石油大學(xué)理學(xué)院,山東青島 266580)

提出一種利用核磁T2譜結(jié)合數(shù)字巖心技術(shù)計算巖石粒度分布的方法。巖石粒度分布近似符合Weibull分布,改變分布參數(shù)可得到不同粒度分布。給出初始Weibull分布后采用過程法構(gòu)建數(shù)字巖心,并用隨機行走法對構(gòu)建的數(shù)字巖心進行核磁數(shù)值模擬,遍歷分布參數(shù)區(qū)間使模擬得到T2譜與試驗T2譜逼近,誤差最小時獲得的Weibull分布即為巖石粒度分布。最后根據(jù)實際核磁測井資料對井下巖石粒度進行連續(xù)計算,粒度計算結(jié)果與實驗室粒度分析結(jié)果對比良好,驗證了方法準確性。

核磁測井；粒度分布；數(shù)字巖心；核磁模擬；核磁T2譜

粒度分析在沉積學(xué)研究中應(yīng)用廣泛,在評述沉積物分類命名方法理論、區(qū)分沉積環(huán)境、判定物質(zhì)輸運方式、判別水動力條件、評價儲層質(zhì)量等方面具有不可替代的作用。早期的粒度分析技術(shù)主要包括直接測量法、篩析法、水析法、激光粒度分析儀法等[1]。但是,這些粒度分析都是在實驗室條件下獲得的,成本較高。近年來國外開發(fā)了一系列新技術(shù)來獲得粒度參數(shù),如超聲波衰減技術(shù)[2]、高清數(shù)字照片自相關(guān)技術(shù)[3]。超聲波衰減技術(shù)只能求取平均粒度,高清數(shù)字照片自相關(guān)技術(shù)井下應(yīng)用受限。而且鉆井取芯僅限于某些井段,不能全面客觀地對儲層進行評價。張旭[4]采用自然伽馬與巖石顆粒的粒度中值建模,實現(xiàn)了井下巖石粒度中值的計算,但這種方法所求的粒度參數(shù)有限。楊寧等[5]提出利用伽馬測井曲線小波變換提取粒度參數(shù),拓展了粒度信息來源,提高了測井資料使用率。筆者提出一種利用核磁T2譜與數(shù)字巖心相結(jié)合連續(xù)計算井下巖石粒度分布的方法。

1 方法原理

巖石的粒度分布可以通過Weibull分布來描述,不同的Weibull分布參數(shù)對應(yīng)不同的粒度分布。巖石的核磁表面弛豫發(fā)生在固液接觸面上,即巖石的顆粒表面,因此核磁T2譜信息能反映巖石的粒度分布。通過選取不同的Weibull分布,采用沉積法生成具有不同粒度分布的數(shù)字巖心,模擬數(shù)字巖心T2譜分布,當(dāng)模擬結(jié)果與試驗T2譜對應(yīng)較好時,所選的Weibull分布可以近似為地層巖石的粒度分布。

流程圖如圖1所示:①確定Weibull參數(shù)取值區(qū)間；②按照指定步長改變Weibull分布參數(shù),得到粒度分布曲線；③根據(jù)粒度分布資料利用過程法和數(shù)學(xué)形態(tài)學(xué)算法構(gòu)建不同含水飽和度的三維數(shù)字巖心并利用隨機行走法模擬數(shù)字巖心核磁T2譜分布；④計算模擬T2譜與試驗T2譜誤差,返回②改變Weibull參數(shù)重新進行模擬計算,直到遍歷所有參數(shù)值,比較計算誤差,誤差最小時對應(yīng)的Weibull分布即為巖石顆粒粒度分布。

圖1 方法流程圖Fig.1 Flowchart of method

1.1 利用Weibull分布確定巖石顆粒粒度分布

研究顯示由于不同的沉積和地質(zhì)環(huán)境,在天然巖石中對于粒度分布有確定的模型,如對數(shù)正態(tài)分布、Weibull分布、對數(shù)雙曲線分布等[6-8]。本次研究中,采用Weibull分布構(gòu)建地層巖石的初始粒度分布。

Weibull概率密度函數(shù):

Weibull累積分布函數(shù):

式中,β為形狀因子；γ為尺度因子；rg為巖石顆粒半徑；rg,0為巖石最小顆粒半徑；P為概率分布；Pc為累積概率分布。

對于不同的β、γ值,Weibull分布形態(tài)是不一樣的,如圖2所示。

圖2 Weibull分布形態(tài)與參數(shù)關(guān)系Fig.2 Relationship between Weibull distribution and parameters

由圖2看出:相同β值下,分布寬度變化較小,隨γ值增大峰值減小,位置右移(圖2(a))；相同γ值下,x坐標值變化較小,隨β值增大峰值增大,且分布范圍縮小(圖2(b))。因此,γ可以近似反映顆粒的粒徑信息的變化,β可以近似反映顆粒的分選信息的變化,實現(xiàn)了數(shù)學(xué)參數(shù)與粒度參數(shù)的對應(yīng)。所以,通過選取不同的參數(shù)值可以獲得不同的粒度分布。圖3是選取的初始粒度分布(β=6,γ=4)。

圖3 初始粒度分布Fig.3 Initial grain size distribution

1.2 應(yīng)用粒度分布構(gòu)建數(shù)字巖心

由Weibull分布獲得的巖石初始粒度分布輸入模型中,根據(jù)孔隙度曲線,采用過程法[9]生成不同孔隙的數(shù)字巖心。其實現(xiàn)的過程如下:

(1)確定一個沉積范圍,使顆粒在給定的范圍內(nèi)沉積。如本文沉積范圍為xy平面0～2500 μm范圍內(nèi)向上堆積。

(2)執(zhí)行下落過程。隨機選取一個下落點即(x,y)確定,使之下落,這時關(guān)鍵是確定z坐標,隨機選取半徑r,所有已存在球的半徑均增加r。

(3)記錄該球下落時所有可能遇到的球,最終將落點最高的球作為落點球。如果沒有落點球則直接到達xy平面z=r。

(4)如果存在落點球,這時落點以球坐標系中與z軸的夾角φ為變量進行移動且要求φ≤90°,同時判定該點是否遇到其他球,如果遇到其他球則將首先遇到的球作為接觸球。否則,從移動的最終位置下落,返回步驟(3)。

(5)在確定落點球與接觸球后,該下落點以二者球心連線作為軸線再次旋轉(zhuǎn),并且轉(zhuǎn)動角度不可以超過一定值,同時確定該點是否遇到其他球。如果遇到其他球,則將其作為第三個球即穩(wěn)定球,同時當(dāng)前位置即為穩(wěn)定位置。否則,從轉(zhuǎn)動的最終位置下落回到步驟(3),繼續(xù)執(zhí)行。

(6)在確定落點和三個球的球心坐標后,判定該點與三個球的關(guān)系是否穩(wěn)定。若穩(wěn)定,則將此時下落點位置作為下落球的球心,所有球的半徑恢復(fù)初值。否則,在這三個球中重新確定落點球和接觸球返回步驟(5),重新確定穩(wěn)定球和穩(wěn)定點。

初始顆粒堆積為巖石壓實和成巖模擬奠定了基礎(chǔ)。壓實過程模擬不考慮顆粒形變和重組效應(yīng),只是將所有顆粒的垂向坐標向模型中心xy平面移動來模擬壓實過程,公式如下:

z=0.5λ(zmax-zmin)+z0(1-λ+ξ).(4)式中,z為新的垂向坐標；z0為初始位置；λ為壓實因子；ξ為模擬顆粒排列的變量；zmax和zmin分別為所有沉積顆粒的最大和最小垂向坐標。圖4(a)是沉積法模擬顆粒堆積成果圖,圖4(b)是對沉積結(jié)果離散化后的成果圖,其中淺灰色表示孔隙,黑色表示骨架。

圖4 過程法構(gòu)建數(shù)字巖心實例Fig.4 Digital rock constructed by process-based method

在100%含水巖石中,T2譜幅度反映了水的含量。然而在部分飽和水(Sw＜100%)的巖石中,核磁T2譜幅度由油水兩相的含量決定。圖5是采用數(shù)學(xué)形態(tài)學(xué)算法[10]構(gòu)建了含水飽和度分別為100%、75%、45%的數(shù)字巖心,其中黑色代表骨架,深灰色代表油、淺灰色代表水。

1.3 數(shù)字巖心NMR弛豫響應(yīng)模擬

儲層巖石通常存在一個孔隙尺寸分布,并且常常含有多種流體成分。因此自旋回波串不是單個T2值的衰減,總的弛豫率為這些弛豫率的疊加[11],

其中,Ai為弛豫時間T2i組分所占的比例,它對應(yīng)于巖石內(nèi)在的比表面(S/V)或孔喉的分布。

核磁共振自旋弛豫的模擬采用隨機行走方法(random walk method)模擬。該算法模擬擴散粒子(稱為walkers或是random walkers)的布朗運動過程,油相和水相核磁模擬過程類似,本文以水相為例簡述模擬過程如下[12-13]:

首先,粒子隨機安放在孔隙空間中。在每個時間步中,粒子從它們的初始位置移動到一個新的位置,該位置處于以初始點為球心半徑為ε(μm)的球表面上。時間步長Δt(ms)定義為

其中D是流體體積擴散系數(shù),μm2/ms。新的位置是

角度θ和φ是隨機選取的。因為粒子在各個方向上行走的概率相等,所以,上面描述的方法僅對于無磁場梯度的系統(tǒng)有效。例如應(yīng)用磁場B0是各向同性的并且孔隙體系中沒有引起磁場梯度的黏土存在。

圖5 不同含水飽和度的數(shù)字巖心Fig.5 Digital rock model with different water saturation

其中,p(t,r)是分布函數(shù),它與隨機粒子移動一個距離r所用的時間t有關(guān)。

然而,對于一個真實的孔隙介質(zhì),構(gòu)建孔隙空間的最大虛構(gòu)球的計算量很大,尤其是每一個時間對于每個粒子都要考慮。在研究中,通常使用半固定步長的方法。本文選擇兩個步長,一個小(ε),另一個步長大,并且為體素尺度的等級。如果一個孔隙點周圍26個點均為孔隙,這是采用大步長同時根據(jù)式(8)進行時間的更新。否則,采用小步長使用式(6)進行時間的更新。

如果粒子進入骨架,它有一定的概率會被“殺死”而消失。消失概率γ與表面弛豫強度ρ有關(guān),使用一個小的固定ε對于大系統(tǒng)的模擬可能會造成計算不夠充分。Zheng和Chiew在每個時間步中使用一個最優(yōu)步長[14]。在每個時間步中,要求最大虛構(gòu)球的球心在粒子初始位置,并且不可以覆蓋任何固體表面。接下來,粒子移動到虛構(gòu)球表面上的任何隨機位置。通過如下初始時間分布進行時間的更新:

如果粒子存活,它的位置不再改變,同時根據(jù)式(6)進行時間的更新。通過對大量的粒子重復(fù)這個過程,可以得到隨機行走粒子的生命分布曲線。這時,這條生命分布曲線可以用來計算磁化強度幅度隨時間的變化曲線。通過式(5)的多指數(shù)反演可以得到T2分布。圖6為試驗T2譜與模擬T2譜結(jié)果對比。圖7為不同含水飽和度數(shù)字巖心模擬T2譜結(jié)果。

圖6 模擬T2譜與試驗T2譜對比Fig.6 Comparison between NMR simulation result and experimental result

圖7 不同含水飽和度核磁模擬結(jié)果Fig.7 NMR simulation result for different water saturation

1.4 參數(shù)最優(yōu)化選擇

不同的Weibull分布參數(shù)對應(yīng)不同的粒度分布,經(jīng)過篩選確定了形狀因子β和尺度因子γ的取值區(qū)間(表1)。以指定步長遍歷形狀因子和尺度因子的值,利用式(10)計算模擬T2譜與試驗T2譜誤差,如果誤差小于10%,則終止遍歷,認為此時確定的分布即為巖石粒度分布。倘若區(qū)間遍歷完成后,計算誤差都大于10%,則采用多組計算整體尋優(yōu)的思想,比較計算誤差,誤差最小時對應(yīng)的參數(shù)選為最優(yōu)化參數(shù),由它們所確定的分布即為巖石粒度分布。

表1 Weibull分布參數(shù)分布區(qū)間Table 1 Parameters? region of Weibull distribution

其中,N是總的布點數(shù),fi,s是點i處模擬T2譜值,fi,m是點i處實測T2譜分布值。

2 處理實例

為了檢驗本文方法的準確性,將本方法計算的巖石粒度結(jié)果與粒度試驗資料作對比,并根據(jù)某井核磁測井資料對井下巖石粒度按照測井采樣間隔進行連續(xù)計算。圖8是采用本文方法計算的粒度分布與試驗測量粒度分布結(jié)果對比圖,對比結(jié)果表明:隨著含水飽和度的降低,粒度分布曲線左移,但是變化較小,計算的粒度分布與實驗結(jié)果相似,可以采用本文方法計算巖石粒度分布。

圖9是采用本文方法連續(xù)計算井下粒度的成果圖,其中第一道是孔隙分析,第二道是核磁T2譜,第三道是計算的粒度概率曲線,第四道是計算的粒度累積概率曲線。右側(cè)是3285.2m,3331 m鉆井取芯粒度試驗測量結(jié)果與計算結(jié)果的對比圖,Weibull分布參數(shù)分別為β=7,γ=4.5與β=10,γ=6.2。試驗結(jié)果與計算結(jié)果有很好的相似性,說明通過核磁T2譜與數(shù)字巖心技術(shù)連續(xù)計算井下巖石粒度分布的方法可行,能夠用來對井下地層粒度參數(shù)進行定性和定量評價。

圖8 粒度計算結(jié)果與含水飽和度關(guān)系Fig.8 Relationship between water saturation and GSD result

圖9 A井粒度計算結(jié)果Fig.9 Grain size distribution calculation result of Well A

3 結(jié) 論

(1)通過核磁T2譜與數(shù)字巖心技術(shù)結(jié)合可以連續(xù)計算井下巖石粒度分布,降低了對鉆井取芯和試驗的依賴,計算結(jié)果與試驗結(jié)果相比具有很好的相似性,驗證了方法的準確性。

(2)通過核磁T2譜與數(shù)字巖心技術(shù)計算巖石粒度分布的方法可行。將該方法應(yīng)用于實際井資料,可以用來對井下地層粒度參數(shù)進行定性和定量評價。

(3)該方法適應(yīng)于砂巖儲層,構(gòu)建的巖石模型只考慮了巖石形成過程中的沉積作用與壓實作用,并未考慮成巖過程。

[1] 肖晨曦,李志忠.粒度分析及其在沉積學(xué)中應(yīng)用研究[J].新疆師范大學(xué)學(xué)報:自然科學(xué)版,2006,25(3): 118-123.

XIAO Chen-xi,LI Zhi-zhong.The research summary of grain size analysis and its application in the sedimentation [J].Journal of Xinjiang Normal University(Natural Sciences Edition),2006,25(3):118-123.

[2] SARPUN I H,SELAMI K K.Mean grain size determination in marbles by ultrasonic velocity techniques[J]. NDT&E International,2005,38(1):21-22.

[3] PATRICK L B,DAVID M R.Field test comparison of an autocorrelation technique for determining grain size using a digital beachball camera versus traditional methods[J]. Sedimentary Geology,2007,201(1/2):180-195.

[4] 張旭.復(fù)雜碎屑巖儲層測井參數(shù)計算方法研究[J].天然氣勘探與開發(fā),2003,26(3):17-23.

ZHANG Xu.Study on methods of computing logging parameters in complex clastic reservoir[J].Natural Gas Exploration and Development,2003,26(3):17-23.

[5] 楊寧,王貴文,賴錦,等.應(yīng)用伽馬測井曲線小波變換計算粒度參數(shù)[J].現(xiàn)代地質(zhì),2012,26(4):778-783.

YANG Ning,WANG Gui-wen,LAI Jin,et al.Application of Gamma curves wavelet transform to calculate grain size parameters[J].Geoscience,2012,26(4):778-783.

[6] BARNDORFF N O,DALSGAARD K,HALGREEN C, et al.Variations in particle size distribution over a small dune[J].Sedimentology,1982,29(1):53-65.

[7] CHRISTIANSEN C,BLAESILD F,DALSGAARD K. Re-interpreting‘segmented’grain size curves[J].Geological Magazine,1984,121(1):47-51.

[8] KONDOLF G M,ADHIKARI A,WEIBULL Vs.Lognormal distributions for fluvial gravels[J].Journal of Sedimentary Research,2000,70(3):456-460.

[9] REN P E,BAKKE S.Process based reconstruction of sandstones and predictions of transport properties[J]. Transport in Porous Media,2002,46(2/3):311-343.

[10] LIU Xue-feng,SUN Jian-meng,WANG Hai-tao.Numerical simulation of rock electrical properties based on digital cores[J].Applied Geophysics,2009,6(1):1-7.

[11] 肖立志.核磁共振成像測井與巖石核磁共振及其應(yīng)用[M].北京:科學(xué)出版社,1998.

[12] ARNS C H,SHEPPARD A P,SOK R M,et al.NMR petrophysical predictions on digitized core images[J]. Petrophysics,2007,48(3):202-221.

[13] HIDAJAT I,SINGH M,COOPER J,et al.Permeability of porous media from simulated NMR response[J]. Transport in Porous Media,2002,48(2):225-247.

[14] ZHENG L H,CHIEW Y C.Computer simulation of diffusion-controlled reactions in dispersions of spherical sinks[J].Journal of Chemical Physics,1989,90(1): 322-327.

(編輯 修榮榮)

Calculation of grain size distribution using NMR T2spectrum and digital rock technology

SUN Jian-meng1,ZHAO Jian-peng1,YAN Wei-chao1,LIU Xue-feng2

(1.School of Geosciences in China University of Petroleum,Qingdao 266580,China；
2.School of Science in China University of Petroleum,Qingdao 266580,China)

A new method calculating grain size distribution(GSD)was proposed by using nuclear magnetic resonance (NMR)T2spectrum and digital rock technology.The curve of GSD is approximately consistent with Weibull distribution,so different GSD can be obtained by changing the parameters of Weibull distribution.According to the initial Weibull distribution,digital rock was constructed by process-based method,and then NMR simulation was conducted on this model by random walk method.To go through the parameter region of Weibull distribution makes the simulated T2spectrum approach the measured one.The Weibull distribution determined by these parameters was the GSD of rock.The method was used to calculate the GSD of underground rock continuously based on NMR logging.The simulated results agree well with the experimental results and validate the accuracy of this method.

nuclear magnetic logging；grain size distribution；digital rock；nuclear magnetic resonance simulation；nuclear magnetic resonance T2spectrum

TE 19

A

1673-5005(2013)03-0057-06

10.3969/j.issn.1673-5005.2013.03.009

2013-02-25

國家科技重大專項(2011ZX05006-002)

孫建孟(1964-),男,教授,博士生導(dǎo)師,主要從事地球物理測井、巖石物理實驗、非常規(guī)油氣測井評價等方面的研究。E-mail: sunjm@upc.edu.cn。