馬百征,許艷爭,黎小偉,龍利平,范久宵

(中國石油化工股份有限公司華北分公司勘探開發(fā)研究院,河南 鄭州 450006)

0 引 言

測井曲線形態(tài)存在的波動變化,可以用信號分析方法對曲線波動變化趨勢分析進(jìn)行巖性識別和地層層序劃分。陳鋼花等[1]用小波時頻分析測井曲線的頻率能量變化進(jìn)行巖性識別和地層層序的劃分,付文釗[2]、李江濤等[3]用小波變換和希爾伯特-黃變換對測井曲線值的突變變化分析進(jìn)行巖性識別和地層層序劃分?？紤]到小波變換、希爾伯特-黃變換等側(cè)重分析樣點(diǎn)數(shù)值的趨勢變化,實(shí)際地層巖性不同是因礦物元素含量、物性結(jié)構(gòu)比例不同,是其樣點(diǎn)值統(tǒng)計的概率分布不同,可以從分析測井?dāng)?shù)值概率分布的角度進(jìn)行巖性識別。

根據(jù)趙志飛等[4]、畢林銳等[5]的著述,可以用正態(tài)分布對地層的巖性、物性分布規(guī)律進(jìn)行表述,某種巖性的測井?dāng)?shù)值在其值域內(nèi)的分布規(guī)律可以認(rèn)為遵循正態(tài)分布。

巖性特征值可以以高斯概率的方式表述,以其巖性值的高斯曲線特征形態(tài)進(jìn)行巖性識別。直觀上可以對于1條測井曲線進(jìn)行逐點(diǎn)逐段的高斯擬合并求出高斯曲線參數(shù)進(jìn)行巖性識別,但該方法分段過小會受到樣點(diǎn)趨勢變化的影響并且計算量會很大。根據(jù)巖性樣點(diǎn)概率統(tǒng)計與樣點(diǎn)數(shù)值趨勢和數(shù)量多少無關(guān)的特性,利用測井曲線整體曲線包絡(luò)是所有局部曲線包絡(luò)疊加的特點(diǎn),可以先計算出某井全井段的巖性包絡(luò),由于全井段的巖性包絡(luò)為各井段不同巖性包絡(luò)的疊加,其擬合出的高斯曲線必定與標(biāo)準(zhǔn)的單一巖性高斯曲線不符,但可以對全井進(jìn)行分段計算包絡(luò)、擬合高斯曲線,當(dāng)分段擬合的高斯曲線符合某種巖性的高斯曲線特征時可以判斷出該井段的巖性。求取某段巖性時要先求取包含該段的上一井段曲線的包絡(luò),計算起始需要給出全井段的曲線包絡(luò),該算法過程可以認(rèn)為是一個遞歸計算過程。

1 實(shí)現(xiàn)方法

1.1 測井曲線直方差包絡(luò)的求取

通常,一段巖性概率值測井曲線的直方差是一段不連續(xù)的折線。以某工區(qū)A井自然伽馬測井曲線為例,圖1中紅色折線為測井曲線的部分井段的直方差曲線,直方差數(shù)值只是1組測量數(shù)據(jù)在有限N次統(tǒng)計的結(jié)果,考慮到隨著N的增加不連續(xù)點(diǎn)之間概率存在的可能性也在增加,需要求出1組光滑的曲線包絡(luò)表示真實(shí)的概率。對于概率曲線,其曲線的極值就是信號的包絡(luò)值。首先找出該曲線的極值點(diǎn),由于極值點(diǎn)存在不等距的情況,需要將其進(jìn)行光滑插值[6-7]得到其包絡(luò)曲線(見圖1藍(lán)色曲線)。

圖1 A井全井段自然伽馬直方差包絡(luò)曲線

圖1中,所有曲線直方差數(shù)值均在其曲線包絡(luò)之內(nèi),包絡(luò)曲線緊密覆蓋直方差曲線,包絡(luò)曲線形態(tài)與直方差曲線形態(tài)一致。直方差的一階包絡(luò)曲線可以認(rèn)為是測井曲線值概率在其值域的最優(yōu)表示。

1.2 測井曲線包絡(luò)的高斯擬合

對于一段井深的測井曲線包絡(luò),其本身形態(tài)的復(fù)雜性和概率分布的偶然性很難僅從其形態(tài)進(jìn)行巖性判定。圖1曲線包絡(luò)中的左側(cè)曲線峰為砂巖峰,右側(cè)為泥巖峰,單從圖1整體的曲線包絡(luò)很難判斷出該段巖性是砂巖還是泥巖、粉砂質(zhì)泥巖、泥質(zhì)砂巖等。根據(jù)測井?dāng)?shù)值遵循的高斯分布規(guī)律,可以對其曲線包絡(luò)進(jìn)行高斯擬合,根據(jù)擬合出的高斯曲線參數(shù)進(jìn)行巖性判斷。標(biāo)準(zhǔn)的高斯分布的概率密度為

(1)

式中,μ為期望值;σ為標(biāo)準(zhǔn)方差。

根據(jù)討論的實(shí)際情況,對于一段自然伽馬測井?dāng)?shù)據(jù)可以將其視為1組離散數(shù)據(jù)(xi,yi)(i=1,2,3,…),其中x為井深,y為自然伽馬值。該數(shù)據(jù)符合高斯分布,可表示為

(2)

(3)

式中,Y和F、S可以理解為高斯曲線的峰值、峰值位置、半寬信息。在本文中峰值位置F、半寬信息S為高斯曲線最重要參數(shù),在后述中將用峰值F表示巖性值、半寬信息S表示巖性的單一性。

對于高斯擬合曲線的參數(shù)F、S值求取常見的如迭代算法[8],但迭代算法需要事先給定模型且計算量很大,如果給定初始模型不合適,難求得正確的結(jié)果。考慮到初始模型給定的復(fù)雜性和迭代的計算量以及對于F、S的求取可采用的算法,對式(3)等式兩邊取自然對數(shù)得到等式

(4)

將式(4)整理

(5)

(6)

簡寫成

A=XB

(7)

式中,A、X矩陣為實(shí)際觀測的已知數(shù)據(jù)矩陣,根據(jù)廣義最小二乘原理可求出最終的矩陣

B=(XTX)-1XTA

(8)

通過矩陣B可以求出峰值位置F、半寬信息S。圖2為某工區(qū)B井部分井段的聲波直方差曲線包絡(luò)(紅色)和擬合出的高斯曲線(藍(lán)色)。

圖2 B井部分井段聲波高斯曲線

1.3 利用高斯曲線參數(shù)進(jìn)行巖性識別方法

圖3 各種巖性的高斯曲線

圖3的(a)、(b)、(c)、(d)分別為各種巖性對應(yīng)的自然伽馬曲線包絡(luò)(紅色)和高斯擬合曲線(藍(lán)色)。圖3(a)為全井段高斯曲線(泥巖砂巖混合F=109,S=32);圖3(b)為砂巖段(含少量泥巖F=100,S=43);圖3(c)為粉砂質(zhì)泥巖(F=110,S=17);圖3(d)為泥質(zhì)砂巖(F=81,S=22)。圖3(a)由于是全井段的自然伽馬,數(shù)值統(tǒng)計井段里含有泥巖和砂巖,曲線包絡(luò)上會出現(xiàn)2個明顯的峰,分別是左側(cè)砂巖峰和右側(cè)泥巖峰;圖3(b)段由于是大段砂巖含少量泥巖,曲線包絡(luò)有明顯的砂巖峰和低幅的泥巖曲線顯示;圖3(c)和圖3(d)均為巖性成分較為單一的粉砂質(zhì)泥巖和砂質(zhì)泥巖,它們的曲線包絡(luò)雖然有距離很近的幾個尖峰,但認(rèn)為是1個峰。通過以上曲線包絡(luò)擬合計算結(jié)果可知,當(dāng)統(tǒng)計的巖性段不為單一巖性時,其曲線包絡(luò)會出現(xiàn)各自對應(yīng)的巖性峰,擬合之后的高斯曲線為幾個峰值的擬合結(jié)果,其曲線幅度較寬,半寬值S會較大,峰值F值反映的是混合巖性值,不能反映真實(shí)的巖性[見圖3(a)、圖3(b)];相反,如果巖性較為單一,則曲線較窄,半寬信息S會很小,峰值F可以反映真實(shí)的巖性[見圖3(c)和圖3(d)]。因此,可以通過井段包絡(luò)的高斯擬合曲線參數(shù)半寬信息S和峰值數(shù)值F識別巖性。

2 應(yīng)用情況

2.1 對區(qū)域已知井標(biāo)準(zhǔn)巖性高斯曲線形態(tài)調(diào)查

在對一個工區(qū)通過測井高斯曲線進(jìn)行巖性識別之前,首先要調(diào)查已知井巖性的測井直方差高斯曲線參數(shù)的分布特征,了解已知巖性和測井直方差高斯曲線參數(shù)之間的關(guān)系。圖4(a)為A井卡本圖上的砂巖段(深度250～350 m),圖4(b)為對應(yīng)的自然伽馬曲線包絡(luò)(紅色)和擬合出的高斯曲線(藍(lán)色),其中圖4(b)圖高斯峰值位置F為59,半寬信息S為8.4。圖4(c)為A井卡本圖上的泥巖段,圖4(d)為該段泥巖自然伽馬高斯曲線,峰值F為125,半寬信息S為18.4。用同樣的方法對該工區(qū)一系列已知井進(jìn)行調(diào)查可以得到該區(qū)域的砂巖泥巖的高斯曲線參數(shù)(見表1)。表1系列測試數(shù)據(jù)顯示,該工區(qū)砂巖峰值位置主要分布在50～90,半寬信息在小于20;泥巖的峰值分布范圍100～150,半寬信息小于35。

表1 某區(qū)系列井砂巖、泥巖測井高斯曲線的調(diào)查情況

圖4 A井泥巖、砂巖與其對應(yīng)的自然伽馬高斯曲線

2.2 用高斯擬合遞歸進(jìn)行巖性識別的過程

圖5 遞歸分離包絡(luò)識別巖性過程

遞歸是利用當(dāng)前結(jié)果逐漸向后遞推計算直到得到最終結(jié)果的計算過程[9],相對于窮舉法,計算效率更高,根據(jù)其算法原理可以設(shè)計出很多精巧高效的程序代碼[10]。本文由于整條測井曲線的巖性包絡(luò)為該井所有巖性包絡(luò)的疊加,需要用遞歸的方法將測井曲線包絡(luò)逐步分離并進(jìn)行高斯曲線擬合。首先對1條測井曲線包絡(luò)進(jìn)行高斯擬合出整條曲線的巖性參數(shù),當(dāng)參數(shù)不滿足單一巖性時繼續(xù)對當(dāng)前井段分段算擬合計算,直到分段的曲線巖性參數(shù)確定為單一巖性時終止分段擬合。理論上,無論泥砂巖怎么交錯,通過逐次對井段的測井曲線包絡(luò)的分離擬合,整個井任意井段的巖性都可以識別出來。

以圖5的A井為例進(jìn)行實(shí)際運(yùn)算。首先將全井段自然伽馬曲線包絡(luò)和擬合高斯曲線求出[見圖5(a)],由于全井段包含泥巖和砂巖,巖性并不單一,擬合出的高斯曲線半寬信息S過大,峰值F所在位置不能反映真實(shí)的巖性,需要繼續(xù)分段求取曲線包絡(luò)和高斯曲線,將A井段的巖性等分,分別求出A井的上半段曲線包絡(luò)和擬合高斯曲線圖5(b)和下半段曲線包絡(luò)和擬合高斯曲線圖5(c);顯然從圖5(b)、圖5(c)上,巖性還不是單一巖性,需要繼續(xù)分段求取包絡(luò)和擬合高斯曲線參數(shù)。繼續(xù)對其中的一段進(jìn)行分段[見圖5(d)、圖5(e)、圖5(f)、圖5(g)、圖5(h)、圖5(i)],直到巖性為單一巖性[見圖5(h)]。圖5(i)中,統(tǒng)計的深度正好處于巖性分界面,包絡(luò)會出現(xiàn)極端數(shù)值過大的2個巖性峰,導(dǎo)致擬合參數(shù)出錯,曲線包絡(luò)無法擬合出其分段正確的高斯曲線,需要繼續(xù)對其分段擬合計算。表2為對某區(qū)的一個標(biāo)準(zhǔn)井?dāng)M合計算情況。

表2中的D1、D2分別為分段的起始深度和結(jié)束深度,運(yùn)算過程是首先對全井段進(jìn)行分段,再選擇其中不符合單一巖性的井段。本文以上段為起算點(diǎn)進(jìn)行逐步分段,直至求出最上段的井巖性值,將剩下的井段當(dāng)作一個未知井段重復(fù)以上過程,最后求出全井段巖性值。通過表2數(shù)據(jù)和已知卡奔圖上的巖性對比,解釋結(jié)果與遞歸高斯擬合分析結(jié)果基本一致,認(rèn)為該方法可以應(yīng)用于基于自然伽馬測井?dāng)?shù)值的巖性識別。

表2 某區(qū)系列實(shí)測A井高斯曲線參數(shù)巖性的調(diào)查情況

3 結(jié)束語

(1) 雖然在理論上可以對某巖性段無限制的分段和包絡(luò)分離,但實(shí)際情況當(dāng)井段深度分的過小、采樣點(diǎn)過少,其擬合出的參數(shù)無法表示正確巖性。

(2) 該方法基于巖性測井值概率統(tǒng)計原理進(jìn)行巖性識別,計算過程中有時會出現(xiàn)不同巖性值混合在一起統(tǒng)計成一種二者相似的巖性,比如單純的砂巖和泥巖,當(dāng)比例一致時有時會識別成砂質(zhì)泥巖或泥質(zhì)砂巖,但這樣的情況往往會伴隨遞歸出的相似巖性分段過大,遇到這樣的情況可以將過大的巖性段繼續(xù)分段,以消除這種情況發(fā)生。

(3) 分段計算過程中遞歸是以均等深度分段的方式計算,應(yīng)用時可以根據(jù)已知地層沉降規(guī)律,巖漿侵蝕、河道運(yùn)移規(guī)律調(diào)試分段深度的比例,以便運(yùn)算更接近地層的真實(shí)分布情況,提高運(yùn)算效率。

(4) 本文僅以單一自然伽馬測井曲線為例進(jìn)行巖性判斷,實(shí)際應(yīng)用中可以同時對其他類型的測井曲線進(jìn)行包絡(luò)分析,實(shí)現(xiàn)更高精度的含油砂巖的識別。比如可以利用自然伽馬、電阻率、聲波、密度、中子等在含油砂巖的曲線特征[11],通過它們同井段包絡(luò)之間的關(guān)系進(jìn)行含油砂巖識別。

(5) 高斯擬合算法是基于地層的巖性、物性遵循高斯分布的情況下進(jìn)行的,對于前期數(shù)據(jù)調(diào)查如果不遵循高斯分布,需要進(jìn)行數(shù)據(jù)整理變換[12]再進(jìn)行高斯擬合,或者根據(jù)其數(shù)據(jù)分布類型進(jìn)行其他類型的擬合,比如泊松擬合、二項(xiàng)式擬合等。

(6) 本文的工作思路是基于整體向局部逐漸遞推的分析方法,不側(cè)重于局部巖性的變化描述。對于巖性值的局部趨勢變化可參考陳鋼花[1]、付文釗[2]、李江濤[3]等用小波變換、希爾伯特-黃變換等分析方法;對于巖性識別和地層層序劃分,不同方法可以相互參考,相互驗(yàn)證。

參考文獻(xiàn):

[1] 陳鋼花,余杰,張孝珍. 基于小波時頻分析的測井層序地層劃分方法 [J]. 新疆石油地質(zhì),2007,28(3): 355-357.

[2] 付文釗,余繼峰,楊鋒杰,等. 小波變換與Hilbert-Huang變換應(yīng)用于層序劃分的比較 [J]. 煤炭學(xué)報,2013,38(2): 435-440.

[3] 李江濤,余繼峰,李增學(xué). 基于測井?dāng)?shù)據(jù)小波變換的層序劃分 [J]. 煤田地質(zhì)與勘探,2004,32(2): 48-50.

[4] 趙志飛,閆暉,姚嵐,等. 正態(tài)分布在區(qū)域地球化學(xué)調(diào)查樣品分析質(zhì)量評價中的應(yīng)用探討 [J]. 巖礦測試,2013,32(1): 97-99.

[5] 畢林銳,毛志強(qiáng),肖承文,等. 正態(tài)分布法在油(氣)水層判別上的應(yīng)用 [J]. 石油天然氣學(xué)報,2006,28(3): 77-78.

[6] 李慶揚(yáng),王能. 超易大義數(shù)值分析 [M]. 北京: 清華大學(xué)出版社,2008.

[7] 楊金才,梁川,盛兆順,等. 數(shù)字信號的包絡(luò)分析方法 [J]. 鄭州工業(yè)大學(xué)學(xué)報,2001,22(3): 59-61.

[8] 王培榮,楊藝. 統(tǒng)計平均擬合正態(tài)分布曲線的迭代算法及其應(yīng)用 [J]. 渝州大學(xué)學(xué)報(自然科學(xué)版),2001,18(1): 37-39.

[9] 邢滔滔. 數(shù)理邏輯 [M]. 北京: 北京大學(xué)出版社,2008.

[10] 譚浩強(qiáng). C程序設(shè)計 [M]. 北京: 北京大學(xué)出版社,2005.

[11] 宋延杰,陳科貴,王向公. 地球物理測井 [M]. 北京: 石油工業(yè)出版社,2011.

[12] 林存山. 地球化學(xué)正態(tài)分布悖論 [J]. 物探化探計算技術(shù),1994,16(4): 289-293.