姜虹，范宜仁，陳華

（1.中國石油大學(xué)地球資源與信息學(xué)院，山東青島266555；2.中國石化股份勝利油田分公司地質(zhì)科學(xué)研究院，山東東營257015；3.中國石油大學(xué)CNPC測井重點實驗室，山東青島266555；4.中國石油大學(xué)理學(xué)院，山東青島266555）

應(yīng)用分形理論計算碳酸鹽巖地層膠結(jié)指數(shù)m值

姜虹1，2，范宜仁1，3，陳華3，4

（1.中國石油大學(xué)地球資源與信息學(xué)院，山東青島266555；2.中國石化股份勝利油田分公司地質(zhì)科學(xué)研究院，山東東營257015；3.中國石油大學(xué)CNPC測井重點實驗室，山東青島266555；4.中國石油大學(xué)理學(xué)院，山東青島266555）

碳酸鹽巖儲層具有復(fù)雜的孔隙空間結(jié)構(gòu)與較強(qiáng)非均質(zhì)性，致使膠結(jié)指數(shù)m值出現(xiàn)了規(guī)律性異化。引入分形幾何理論，以巖心分析數(shù)據(jù)和測井?dāng)?shù)據(jù)為依據(jù)，開展復(fù)雜碳酸鹽巖儲層分形特征及應(yīng)用技術(shù)研究。以巖心數(shù)據(jù)作為一個時間序列，通過嵌入維數(shù)與無標(biāo)度區(qū)的自適應(yīng)計算方法，在完善關(guān)聯(lián)維數(shù)算法的同時，建立了巖心分維值與膠結(jié)指數(shù)m的基本關(guān)系；分析測井曲線的分形特性，以分形理論為橋梁，建立了測井曲線分維值與m值之間的函數(shù)關(guān)系，實現(xiàn)了通過測井曲線求取m值的設(shè)想。

測井曲線；分形理論；關(guān)聯(lián)維數(shù)；巖心分析；膠結(jié)指數(shù)；碳酸鹽巖

0 引 言

分形［1］幾何在石油地質(zhì)、巖石物性等方面已經(jīng)得到廣泛應(yīng)用。近期研究表明，大多數(shù)巖石的原生、次生孔隙，以及裂縫和溶洞網(wǎng)絡(luò)，即使其孔隙尺寸變化極大（納米級至幾十米級尺度范圍），但均呈現(xiàn)分形結(jié)構(gòu)的特點。這些分形對象的大小、形狀、位置以及孔隙度、滲透率等其他屬性參數(shù)的統(tǒng)計分布，都將遵從儲層非均質(zhì)性指數(shù)定律，可用分形幾何進(jìn)行描述，同時這些貌似復(fù)雜的現(xiàn)象僅用少量的參數(shù)（如分?jǐn)?shù)維）便可以表征［2］。

分形理論在測井中也有一定的應(yīng)用，分?jǐn)?shù)維與阿爾奇公式的膠結(jié)指數(shù)m值直接相關(guān)，且已有利用測井曲線資料與地層的分形特性求取地層膠結(jié)指數(shù)m值的先例。鑒于巖心實驗數(shù)據(jù)能夠準(zhǔn)確且直觀反映地層的孔滲特征，可以從巖心數(shù)據(jù)出發(fā)，研究巖石孔隙結(jié)構(gòu)的分形特征。本文分別從巖心數(shù)據(jù)和測井曲線數(shù)據(jù)2個角度，利用分形理論建立分?jǐn)?shù)維與膠結(jié)指數(shù)m值的基本關(guān)系。

1 關(guān)聯(lián)維數(shù)算法

1.1 G－P算法

關(guān)聯(lián)維數(shù)的計算方法有很多種，其中應(yīng)用比較廣泛的是G－P算法。G－P算法基本過程［3］：① 相空間重構(gòu)：將一個指定的時間或空間數(shù)據(jù)系列構(gòu)造成一個d維的相空間［4］（d稱為嵌入維數(shù)）；② 求出以上向量當(dāng)中所有點對的距離γij；給出一個數(shù)r＞0，將所有點對的距離γij與r相比較，求出小于r的γij數(shù)目；③求關(guān)聯(lián)積分函數(shù)C（r）：所有小于r的點對距離γij的數(shù)目在所有點對數(shù)中的比例。

若分形存在，r取充分小值，在某一區(qū)間內(nèi)關(guān)聯(lián)積分函數(shù)會逼近式（1）。C（r）＝rD，其中D為關(guān)聯(lián)維數(shù)。所以，將關(guān)聯(lián)維數(shù)定義為

在上述算法過程中，核心問題是嵌入維數(shù)d和無標(biāo)度區(qū)的確定。如圖1所示，不同的嵌入維數(shù)所對應(yīng)的ln r－ln C（r）的雙對數(shù)曲線不同，導(dǎo)致關(guān)聯(lián)維數(shù)計算不穩(wěn)定、不準(zhǔn)確。而無標(biāo)度區(qū)的確定通常是人工給出，缺乏客觀性。有必要研究適合巖心分析數(shù)據(jù)的嵌入維數(shù)d和無標(biāo)度區(qū)的自動計算方法。

圖1 嵌入維數(shù)1：20對應(yīng)的ln r－ln C（r）的雙對數(shù)曲線圖

1.2 嵌入維數(shù)d的確定

一般情況下，在給定時間序列進(jìn)行相空間重構(gòu)時，其嵌入維數(shù)的確定根據(jù)Takens定理［5］得到，認(rèn)為當(dāng)嵌入維數(shù)d≥2D＋1（D為原動力系統(tǒng)吸引子維數(shù)）時，重構(gòu)相空間和原動力學(xué)系統(tǒng)的相空間拓?fù)涞葍r。但這一結(jié)論的前提是數(shù)據(jù)量無限長且無噪聲干擾，顯然，實際測量的巖心分析數(shù)據(jù)無法滿足以上要求，而且計算的嵌入維數(shù)未必是最小嵌入維數(shù)，這樣容易增加計算的復(fù)雜度，影響計算的準(zhǔn)確性。因此選擇一個合適的嵌入維數(shù)，進(jìn)行有效的相空間重構(gòu)十分關(guān)鍵。

使用Cao方法確定嵌入維數(shù)，定義一個E（d）［6］

式中，d為嵌入維數(shù)；N為數(shù)據(jù)列長度；a（i，d）＝‖yi（d＋1）－yn（i，d）（d＋1）‖／‖yi（d）－yn（i，d）（d）‖，yi（d）為重構(gòu)相空間中第i個向量，yi（d＋1）是以d＋1為嵌入維數(shù)重構(gòu)的相空間中的第i個相量，n（i，d）（1≤n（i，d）≤N－d）是一個整數(shù)，yn（i，d）（d）是在重構(gòu)的d維相空間中與yi（d）最鄰近的點，如果yn（i，d）（d）＝y(tǒng)i（d），將用第2鄰近的點代替它。

隨著d不斷增大，E1（d）將會趨近于一個穩(wěn)定值（程序中設(shè)定了一個下限值δ，當(dāng)ΔE1（d）＜δ時就認(rèn)為E1（d）趨近于一個定值），當(dāng)E1（d）不隨d的增大而變化時，就將此時的d定義為重構(gòu)相空間的嵌入維數(shù)（見圖2，嵌入維數(shù)d＝10）。圖3為確定嵌入維數(shù)為10對應(yīng)的ln r－ln C（r）雙對數(shù)曲線圖。在確定嵌入維數(shù)之后，無標(biāo)度區(qū)的確定通常是人工給出的，缺乏客觀性，因而有必要研究適合巖心分析數(shù)據(jù)的無標(biāo)度區(qū)的自動計算方法。

圖2 嵌入維數(shù)d－E1（d）關(guān)系圖

圖3 嵌入維數(shù)為10的ln r－ln C（r）曲線圖

1.3 無標(biāo)度區(qū)的確定

一般研究的分形系統(tǒng)并非嚴(yán)格意義上的自相似，嚴(yán)格的自相似僅存在于一定的尺度變化范圍內(nèi)，一旦超出了這個尺度變化范圍，其自相似性就不存在，這個尺度的變化范圍就是分形無標(biāo)度區(qū)。根據(jù)分形理論的標(biāo)度不變性，當(dāng)r在分形的無標(biāo)度區(qū)內(nèi)變換時，函數(shù)C（r）與r之間呈指數(shù)關(guān)系，即ln r與ln C（r）之間呈線性關(guān)系。

①根據(jù)G－P算法，得到ln r－ln C（r）的雙對數(shù)點圖，將所有的點擬合成較平滑的n次曲線，并對該曲線上的每個點分別求一階導(dǎo)數(shù)和二階導(dǎo)數(shù)（見圖4）。②由于理想的線性區(qū)域的二階導(dǎo)數(shù)為0，因此為了尋找到無標(biāo)度區(qū)域，定義了一個門限值ε，將曲線上二階導(dǎo)數(shù)值小于ε的每個區(qū)間都找出來。③分別將每個區(qū)間中的點進(jìn)行最小二乘法擬合，并求出擬合度以及擬合后直線的斜率。④ 比較各區(qū)間的點擬合成直線后的擬合度，擬合程度最高的區(qū)間即為無標(biāo)度區(qū)，該區(qū)間直線斜率為關(guān)聯(lián)維數(shù)（見圖5）。

圖4 ln r－ln C（r）擬合曲線及一、二階導(dǎo)數(shù)圖

圖5 無標(biāo)度區(qū)擬合度99.16%

2 豪斯道夫維數(shù)算法

具有分形特征的圖形其覆蓋尺度與測度之間存在指數(shù)關(guān)系，測井曲線為采樣間隔相等的曲線，定義覆蓋曲線的尺度ε與曲線測度μ（k，ε）的關(guān)系式為［7－9］

式中，Δh為測井曲線采樣點的采樣間隔；x（i）為深度i·Δh的測井值；ε＝M·Δh為覆蓋測井曲線的尺度；μ（k，ε）為點k處對應(yīng)尺度ε的測度，當(dāng)覆蓋尺度ε縮?。ɑ驍U(kuò)大）a倍，測度μ（k，ε）將縮?。ɑ驍U(kuò)大）b倍，測井曲線的豪斯道夫維數(shù)Df可表示為［10－11］

當(dāng)豪斯道夫維數(shù)Df大于拓?fù)渚S數(shù)時，可以認(rèn)為曲線具備分形特征。測井曲線具備分形特征的條件即分形維數(shù)Df應(yīng)在1～2之間變化。

3 相似維數(shù)算法

如果一分形集合S可以劃分為N個同等大小的子集，每一子集為原集合放大r倍，則集合S的相似維數(shù)dS定義為［3］

相似維數(shù)值范圍在一般曲線維數(shù)值1和平面維數(shù)值2之間。

4 將關(guān)聯(lián)維數(shù)算法應(yīng)用于巖心數(shù)據(jù)求取膠結(jié)指數(shù)m值

選取某碳酸鹽巖區(qū)塊的5口井進(jìn)行研究，以這5口井取心的巖心分析數(shù)據(jù)為研究對象，以巖心實驗測量的孔隙度為原始數(shù)據(jù)進(jìn)行處理，依據(jù)上述關(guān)聯(lián)維數(shù)的計算方法，分別對每一口井取心的孔隙度進(jìn)行了關(guān)聯(lián)維數(shù)的計算分析，結(jié)果見表1。

表1 某區(qū)塊5口井膠結(jié)指數(shù)m與關(guān)聯(lián)維數(shù)D數(shù)據(jù)表

以此建立的模型方程m＝2.4251－0.2119D，擬合程度為95.0%。將×6井的數(shù)據(jù)代入建立的模型進(jìn)行驗證，×6井的取心孔隙度的關(guān)聯(lián)維數(shù)值為1.396，其膠結(jié)指數(shù)m值為2.149。將×6井的取心孔隙度的關(guān)聯(lián)維數(shù)值為1.396代入模型計算得到的m0＝2.16，誤差分析（m－m0）／m≈0.005。在誤差允許范圍內(nèi)，在該區(qū)塊內(nèi)，上述模型比較適用。

將以上模型與孔隙度－膠結(jié)指數(shù)的交會圖對比（見表2），得出其擬合程度僅為86.3%，對比后可以明顯看出膠結(jié)指數(shù)m與孔隙度φ之間沒有明顯的線性相關(guān)性。將計算得到的孔隙度的關(guān)聯(lián)維數(shù)D與膠結(jié)指數(shù)m交會，能得出兩者有一定線性相關(guān)性。

表2 某區(qū)塊5口井的膠結(jié)指數(shù)與孔隙度數(shù)據(jù)表

以上是對該區(qū)塊各井取心的孔隙度數(shù)據(jù)進(jìn)行的關(guān)聯(lián)維數(shù)計算，得到孔隙度的關(guān)聯(lián)維數(shù)與膠結(jié)指數(shù)m值之間有一定的線性相關(guān)性。說明碳酸鹽巖地層的孔隙應(yīng)該具有一定的分形特征，膠結(jié)指數(shù)m值變化范圍大，但仍有規(guī)律可循。

5 分形理論應(yīng)用于測井曲線計算m值

利用計算得到的測井曲線分維數(shù)就可以間接反映地層信息，如地層的孔隙結(jié)構(gòu)、非均質(zhì)性等地層信息都能夠比較準(zhǔn)確求出。

5.1 測井曲線的分形特征分析

圖6至圖9分別為聲波、電阻率、自然伽馬及井徑測井曲線的分形特征圖。從圖6至圖9可以看出，聲波曲線（Df＝1.03）、電阻率曲線（Df＝1.13）、自然伽馬曲線（Df＝1.04）以及井徑曲線（Df＝1.01）的豪斯道夫維數(shù)都大于其拓?fù)渚S數(shù)1，而且擬合程度都在97%以上，因此以上測井曲線具備分形特征，可以用來進(jìn)行分形計算。

圖6 聲波測井（AC）曲線分形特征

5.2 相似維數(shù)法求m值

在一段儲層中，用邊長為L1的正方形網(wǎng)格覆蓋電阻率測井曲線，邊長為L2的正方形網(wǎng)格覆蓋聲波測井曲線，可以得到2種測井曲線的相似維數(shù)d（Rt）和d（lnφ）。用阿爾奇公式進(jìn)行m值計算時，如果該段儲層的流體性質(zhì)以及巖石的巖性物性變化不大，可以推導(dǎo)得到利用測井曲線分?jǐn)?shù)維計算膠結(jié)指數(shù)m值（令電阻率測井曲線的覆蓋尺度L1與聲波測井曲線的覆蓋尺度L2都等于L）的表達(dá)式［12］

利用上述方法，對某碳酸鹽巖區(qū)塊××2井進(jìn)行膠結(jié)指數(shù)m值計算，取××2井4 921.4～5 203.36m井段，按照儲層的巖性、物性以及流體性質(zhì)特征，將該井段分為11個層段，取L為10－4分別對聲波測井曲線和深側(cè)向測井曲線進(jìn)行覆蓋，得到lnφ和ln Rt，進(jìn)而進(jìn)算出膠結(jié)指數(shù)m值，并將計算值與巖電實驗測量結(jié)果進(jìn)行對比（見表3）。

圖7 電阻率測井（Rt）曲線分形特征

圖8 自然伽馬測井（GR）曲線分形特征

圖9 井徑測井（CAL）曲線分形特征

將通過分形計算得到的膠結(jié)指數(shù)m值與取心實驗測量結(jié)果對比可見，平均誤差在8.0%以下，在誤差允許范圍內(nèi)，計算結(jié)果能夠比較準(zhǔn)確地反映巖石的孔隙結(jié)構(gòu)特征，因此利用測井曲線分維數(shù)求膠結(jié)指數(shù)m值可行。常規(guī)采用巖石物理實驗方法求取膠結(jié)指數(shù)m最直接而且準(zhǔn)確，但是取心和巖石物理實驗過程復(fù)雜而且耗費(fèi)人力物力，通過測井曲線計算得到m值，是比較有意義的。

表3 利用測井曲線分維數(shù)計算膠結(jié)指數(shù)m值與實驗測量值對比表

5.3 關(guān)聯(lián)維數(shù)法求m值

測井曲線采樣點也是一個時間序列，將關(guān)聯(lián)維數(shù)應(yīng)用于測井曲線求取m值方法同樣適用。為建立巖心m值與測井曲線關(guān)聯(lián)維數(shù)值的一一對應(yīng)關(guān)系，需將巖心歸位后，求取歸位深度處各相關(guān)測井曲線的關(guān)聯(lián)維數(shù)值。將每個測井?dāng)?shù)據(jù)序列求出的關(guān)聯(lián)維數(shù)值，作為該數(shù)據(jù)列中點數(shù)據(jù)深度處的分維值，基本能夠反映該深度處的測井曲線分維特征，并且做到了分維值與深度的一一對應(yīng)。

常規(guī)測井曲線當(dāng)中，與孔隙結(jié)構(gòu)指數(shù)m值密切相關(guān)的主要是三孔隙度曲線（聲波、中子和密度），其分形維數(shù)可表示為dAC、dCNL、dDEN，由疊加和匹配原理可以得到一個與地層孔隙結(jié)構(gòu)相關(guān)的分形維數(shù)值dm

［13］［見式（8）］，dm綜合反映了地層信息，與膠結(jié)指數(shù)m之間應(yīng)該有相關(guān)性。

將巖心歸位后，首先計算歸位深度處的dm值，將dm值與巖心分析得到的膠結(jié)指數(shù)m之間進(jìn)行比較（見表4），分別取某碳酸鹽巖區(qū)塊某組上、下段巖心與測井曲線分維值進(jìn)行對比分析，發(fā)現(xiàn)dm與m之間具有較好的線性關(guān)系，見圖10與圖11，擬合關(guān)系式見式（9）和式（10）。

表4 利用測井曲線關(guān)聯(lián)維數(shù)求取膠結(jié)指數(shù)m值

圖10 某組上段m與dm關(guān)系圖

某組上段m與dm的關(guān)系式

某組下段m與dm的關(guān)系式

圖11 某組下段m與dm關(guān)系圖

式（9）、式（10）的分形幾何規(guī)律與儲層的地質(zhì)特征是相符的。分維數(shù)本身反映的是樣本空間的變化劇烈程度，dm綜合反映三孔隙度曲線的變化劇烈程度，dm越大則地層結(jié)構(gòu)越復(fù)雜，即dm與地層的非均質(zhì)性成正比。儲層級別越高，其孔隙的連通性越好，非均質(zhì)性相對較弱，dm越?。环粗?，則dm越大。通過對某碳酸鹽巖區(qū)塊巖心實驗分析得到，高級別儲層的m值較大（原因是高級別儲層主要以孔洞型為主），因此dm與m值理論上應(yīng)該成反比關(guān)系，與推導(dǎo)的關(guān)系式相符。為驗證關(guān)系式的正確性，由式（9）、式（10）分別對×6井的3個類儲層進(jìn)行分形計算，得到各類儲層的m值與巖心實驗分析得到的m值較相符，具體見表5。

表5 巖心膠結(jié)指數(shù)m值按儲層分類統(tǒng)計表

膠結(jié)指數(shù)m值與測井曲線的分維值建立良好的線性關(guān)系，尤其針對強(qiáng)非均質(zhì)性的碳酸鹽巖儲層，m值變化大，且變化規(guī)律不易把握，通過引入分形數(shù)學(xué)的方法，降低了m值求取難度，同時證明復(fù)雜的碳酸鹽巖儲層孔隙結(jié)構(gòu)遵循非均質(zhì)性指數(shù)定律，可用分形幾何進(jìn)行描述。

6 結(jié) 論

（1）實現(xiàn)了關(guān)聯(lián)維數(shù)的自動求取方法，尤其解決了適合對巖心分析數(shù)據(jù)進(jìn)行分維數(shù)計算的嵌入維數(shù)以及無標(biāo)度區(qū)的自動計算方法。

（2）將巖心分析孔隙度數(shù)據(jù)轉(zhuǎn)化為時間序列，進(jìn)行了關(guān)聯(lián)維數(shù)的計算，得到孔隙度的分維數(shù)與膠結(jié)指數(shù)m值之間有一定的線性相關(guān)性，說明碳酸鹽巖儲層的孔隙應(yīng)該具有一定的分形特征。

（3）證明測井曲線具有分形特征，并以測井曲線的分維數(shù)為橋梁，建立了測井曲線與m值之間的函數(shù)關(guān)系，實現(xiàn)了通過測井曲線求取m值的設(shè)想。

［1］ Mandelbrot B B.The Fractal Geometry of Nature［M］.New York：W.H.Freeman And Company，1982.

［2］ 牟書令.中國海相油氣勘探理論技術(shù)與實踐［M］.北京：地質(zhì)出版社，2009.

［3］ 王域輝，廖淑華.分形與石油［M］.北京：石油工業(yè)出版社，1994.

［4］ 陳颙.分形與混沌在地球科學(xué)中的應(yīng)用［M］.北京：學(xué)術(shù)期刊出版社，1989.

［5］ Takens M.Detecting Strange Attractors in Fluid Turbulence［C］∥Dynamical Systems and Turbulence，Lecture Notes in Mathematics.Berlin：Springer.1986.

［6］ Liangyue Cao.Practical Method for Determining the Minimum Embedding Dimension of a Scalar Time Series［J］.Physica D，1997，110：43－50.

［7］ Pape H，Clause C，Iffland J.Variation of Permeability with Porosity in Sandstone Diagenesis Interpreted with a Fractal Pore Space Model［J］.Pure and Applied Geophysics，2000，157：603－619.

［8］ Cheng Q M.Multifractal Interpolation［C］∥Lippard S J，Nass A，Sinding L.5th Annual Conference［LAMG99］Proceedings：International Association of Mathematical Geology.Norway：Trondheim.1999，245－250.

［9］ 李慶謀，成秋明.測井曲線分形校正［J］.中國地質(zhì)大學(xué)學(xué)報：地球科學(xué)，2002，27（1）：63.

［10］謝季堅，鄧小炎.現(xiàn)代數(shù)學(xué)方法選講［M］.北京：高等教育出版社，2001.

［11］梁齊瑞，趙世龍，丁忙生.自然伽馬測井曲線的分形特征分析［J］.物探與化探，2006，30（3）：240－243.

［12］劉紅岐，夏宏泉，王擁軍，等.地層膠結(jié)指數(shù)m的分形特征研究［J］.測井技術(shù)，2001，25（1）：24－27.

［13］王天波，董春旭，徐麗萍，等.用分形理論確定m、n值的方法及其應(yīng)用［J］.測井技術(shù)，1998，22（1）：16－19.

On Calculation of Cementation Index min Carbonate Formation Using Fractal Theory

JIANG Hong1，2，F(xiàn)AN Yiren1，3，CHEN Hua3，4
（1.College of Geo－resource and Information，China University of Petroleum，Qingdao，Shandong 266555，China；2.Geoscience Research Institute，Shengli Oilfield CO.，SINOPEC，Dongying，Shandong 257015，China；3.CNPC Key Well Logging Laboratory，China University of Petroleum，Qingdao，Shandong 266555，China；4.College of Science，China University of Petroleum.Qingdao，Shandong 266555，China）

Carbonate reservoir with complex pore structure and strong heterogeneity often results in irregular change of the cementation index m.So，fractal geometry theory is introduced.Based on core analysis and log data，carried out are the researches on fractal characteristics and application technology in complex carbonate reservoirs.The core data are transformed into time series，and by using adaptive method of embedding dimension and scale－free zone，the correlation dimension algorithm is modified and meanwhile the basic relationship between fractal dimension of core and the cementation index mis established.Then analyzed are the fractal characteristics of logging curves and built is a function of the fractal dimension of logging curves and mthrough the fractal theory，and therefore figuring out the value mwith logging curves.

logging curve，fractal theory，correlation dimension，core analysis，cementation index，carbonate

P631.84；TE19

A

2012－01－05 本文編輯 余迎）

1004－1338（2012）03－0250－06

姜虹，女，1984年生，碩士，研究方向為測井理論、方法與技術(shù)。