姜虹，范宜仁，陳華

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

應用分形理論計算碳酸鹽巖地層膠結指數m值

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

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

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

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

0 引 言

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

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

1 關聯(lián)維數算法

1.1 G－P算法

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

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

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

圖1 嵌入維數1：20對應的ln r－ln C（r）的雙對數曲線圖

1.2 嵌入維數d的確定

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

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

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

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

圖2 嵌入維數d－E1（d）關系圖

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

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

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

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

圖4 ln r－ln C（r）擬合曲線及一、二階導數圖

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

2 豪斯道夫維數算法

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

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

當豪斯道夫維數Df大于拓撲維數時，可以認為曲線具備分形特征。測井曲線具備分形特征的條件即分形維數Df應在1～2之間變化。

3 相似維數算法

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

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

4 將關聯(lián)維數算法應用于巖心數據求取膠結指數m值

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

表1 某區(qū)塊5口井膠結指數m與關聯(lián)維數D數據表

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

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

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

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

5 分形理論應用于測井曲線計算m值

利用計算得到的測井曲線分維數就可以間接反映地層信息，如地層的孔隙結構、非均質性等地層信息都能夠比較準確求出。

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

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

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

5.2 相似維數法求m值

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

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

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

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

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

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

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

5.3 關聯(lián)維數法求m值

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

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

［13］［見式（8）］，dm綜合反映了地層信息，與膠結指數m之間應該有相關性。

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

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

圖10 某組上段m與dm關系圖

某組上段m與dm的關系式

某組下段m與dm的關系式

圖11 某組下段m與dm關系圖

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

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

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

6 結 論

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

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

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

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

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

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

［4］ 陳颙.分形與混沌在地球科學中的應用［M］.北京：學術期刊出版社，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］.中國地質大學學報：地球科學，2002，27（1）：63.

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

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

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

［13］王天波，董春旭，徐麗萍，等.用分形理論確定m、n值的方法及其應用［J］.測井技術，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年生，碩士，研究方向為測井理論、方法與技術。