一種基于巖石薄片圖像的粒度分析新方法

袁　瑞，朱　銳，瞿建華，孫玉秋，唐　勇，潘　進

（1.長江大學地球科學學院，武漢430100；2.長江大學信息與數學學院，湖北荊州434023；3.中國石油新疆油田分公司勘探開發研究院，新疆克拉瑪依834000）

一種基于巖石薄片圖像的粒度分析新方法

袁瑞1，2，朱銳1，瞿建華3，孫玉秋2，唐勇3，潘進1

（1.長江大學地球科學學院，武漢430100；2.長江大學信息與數學學院，湖北荊州434023；3.中國石油新疆油田分公司勘探開發研究院，新疆克拉瑪依834000）

鑒于巖石薄片粒度分析的實驗室技術具有較大的局限性，且速度緩慢，提出了一種基于巖石薄片圖像空間自相關系數的粒度分析新方法。該方法通過計算已知顆粒大小的巖石薄片的空間自相關系數，并利用帶約束條件的最小二乘法求解對應的空間自相關系數方程組，從而得到顆粒大小未知的巖石薄片的粒度分布。采用橢圓隨機生成的方法制作了顆粒大小不同的理論模擬巖石薄片，并分析了其空間自相關系數與粒度的關系，得出空間自相關系數隨著偏移距離的增大而減小；當偏移距離一定時，空間自相關系數隨著粒度的增大而增大。采用該方法對粗砂巖中粒徑分別為0.5～1.0 mm，1.0～1.5 mm和1.5～2.0 mm的顆粒進行了粒度分布的定量計算，得出這3種粒徑顆粒所占的百分比分別為55.8%，24.6%和20.2%，該結果與實際值較接近，變化趨勢與實際一致。

空間自相關系數；粒度分析；最小二乘法；模擬巖石薄片

0　引言

巖石薄片粒度分析是沉積學研究的基礎［1-3］。在取得巖石樣品后，研究人員一般通過自己的經驗進行判斷，這樣的結果往往因人而異，缺乏可靠的數據證明；或者通過實驗進行分析和測量，但需要較長的時間；或者是在巖石薄片上手動測量各個顆粒的直徑，再統計薄片的粒度分布，工作繁瑣。筆者提出一種基于巖石薄片圖像空間自相關系數的粒度分析新方法，采用橢圓隨機生成的方法制作大量已知顆粒大小的理論模擬巖石薄片，分析其空間自相關系數與粒度的關系，最后利用帶約束條件的最小二乘法求解空間自相關系數方程組，從而定量計算出粒度分布。

1　圖像的空間自相關系數定義

同一幅圖像中2個矩形區域間的互相關程度表征了圖像的空間自相關系數，也就是說圖像的空間自相關系數由一系列的互相關系數組成。將巖石薄片的圖像經過預處理后形成灰度圖像。根據統計學原理，巖石薄片灰度圖像內大小相同的X和Y區域內的互相關系數［4］為

式中：R（X，Y）≤1；n為巖石薄片圖像中X和Y區域的橫向寬度，pixel；BX（i，j）和BY（i，j）分別為X和Y區域中橫向上第i個、縱向上第j個像素的灰度值；BX和BY分別為X和Y區域中所有像素灰度的平均值。

為了得到空間自相關系數，引入偏移距離k（以像素為單位，k=0，1，…，K），其定義為圖像中像素大小相同的X和Yk區域間的水平方向距離。巖石薄片圖像空間自相關系數的計算原理如圖1所示，其中X區域右邊界與整個圖像的右邊界相距K個像素點，X區域固定不動。圖像自相關系數的計算步驟如下：①取k=0，區域Y0與X區域重合，利用式（1）計算X與Y0的互相關性系數R（X，Y0）；②取k= 1，2，…，將區域Yk逐漸向巖石薄片圖像的右邊界移動，利用式（1）逐步計算出X與Yk區域的互相關系數R（X，Yk）；③當k=K時，利用式（1）計算出X與YK區域的互相關系數R（X，YK），計算結束。最后令R（k）=R（X，Yk）（k=0，1，…，K），即可得到巖石薄片圖像最大偏移距為K時的空間自相關系數。當X與Y區域的偏移距離k=0時，兩者具有最大的互相關性，即R（X，Y0）=1。

圖1　巖石薄片圖像空間自相關系數計算示意圖Fig.1　Schematic diagram of spatial autocorrelation coefficient computation for rock section

2　粒度分布的定量計算方法

利用空間自相關系數法，可根據已知粒度分布的巖石薄片的空間自相關系數，定量計算顆粒大小未知的巖石薄片的粒度分布。取最大偏移距離為K，計算M種粒徑巖石薄片（顆粒大小已知）的空間自相關系數。在偏移距離k處（k=1，2，…，K），空間自相關系數記為a（k，m）（m=1，2，…，M），顆粒大小未知的巖石薄片的空間自相關系數記為bk，設其M種粒徑的百分比為xm（0≤xm≤1），則滿足以下線性方程組［5］：

當K=M時，式（2）可通過矩陣代數來求解。然而，在實際的巖石薄片圖像中，一般K＞M，因此，采用帶約束條件的最小二乘法來求解，即可得到顆粒大小未知的巖石薄片中M種粒徑的粒度分布。

3　粒度分析實例

3.1空間自相關系數與粒度的關系

為了從理論上分析巖石薄片圖像空間自相關系數法在顆粒粒度分析中的可行性，采用橢圓隨機生成的方法［6］，在10 mm×10 mm的范圍內模擬了不同粒徑的巖石薄片圖像（圖2）。圖2中巖石顆粒的灰度值為255（白色），填充物或泥質的灰度值為0（黑色）。

粗砂巖（0.5～2.0 mm）的模擬薄片中共有62個巖石顆粒（橢圓顆粒），粒徑為0.51～1.96 mm，平均為0.9 mm［圖2（a）］；中砂巖（0.25～0.50 mm）的模擬薄片中共有464個巖石顆粒，粒徑為0.25～0.49 mm，平均為0.34 mm［圖2（b）］；細砂巖（0.05～0.25 mm）的模擬薄片中共有2398個巖石顆粒，粒徑為0.05～0.25 mm，平均為0.14 mm［圖2（c）］。

圖2　理論模擬巖石薄片Fig.2　Simulated theoretical rock sections

如圖3所示，取最大偏移距離K為100時，分別計算出粗砂巖、中砂巖和細砂巖模擬薄片圖像的空間自相關系數。從圖3可以看出，巖石薄片的空間自相關系數隨著偏移距離k的增大而減小；當偏移距離k一定時，空間自相關系數隨著粒度的增大而增大［7-10］。

圖3　砂巖模擬薄片圖像空間自相關系數Fig.3　Spatial autocorrelation coefficient of sandstone sections

圖4　不同粒徑粗砂巖的模擬巖石薄片Fig.4　Simulated gritstone sections with different particle size

3.2粒度分布定量計算

選取圖2（a）中粗砂巖的3種粒徑0.5～1.0 mm，1.0～1.5 mm和1.5～2.0 mm進行粒度分布的定量計算，分別模擬了對應顆粒大小的薄片（圖4）。取最大偏移距離K為100，分別計算得到了這3種粒徑的粗砂巖模擬薄片圖像的空間自相關系數（圖5）。根據空間自相關系數，采用帶約束條件的最小二乘法計算得到這3種粒徑的顆粒所占百分比分別為55.8%，24.6%和20.2%（表1）。這3種粒徑顆粒實際所占的百分比分別為61.3%，32.3%和6.5%。計算結果與實際值接近，總體變化趨勢與實際一致。

圖5　粗砂巖模擬薄片圖像空間自相關系數Fig.5　Spatial autocorrelation coefficient of gritstone sections

表1　粒度分布實際值與計算值對比Table 1　Comparison of grain size distribution between calculated and actual values

4　結論與認識

（1）提出了一種基于巖石薄片圖像空間自相關系數的粒度分析新方法。計算已知顆粒大小的巖石薄片圖像的空間自相關系數，利用帶約束條件的最小二乘法求解空間自相關系數方程組，從而得到顆粒大小未知的巖石薄片中的粒度分布。

（2）利用橢圓隨機生成的方法模擬了巖石薄片，得出其空間自相關系數隨著偏移距離的增大而減小；當偏移距離一定時，空間自相關系數隨著粒度的增大而增大。

（3）利用這種新方法對一組模擬的巖石薄片進行了粒度分析，計算結果與實際值較接近。因此，在取得高質量的巖石薄片后，可嘗試利用該方法進行粒度分析的定量計算。

（References）：

［1］馬洪濤，蔡明，付國民.塔河油田AT1井區中三疊統中油組湖底扇沉積特征研究［J］.巖性油氣藏，2010，22（3）：53-58. Ma Hongtao，Cai Yue，Fu Guomin.Sedimentary characteristics of sublacustrine fan of Middle Triassic Zhongyou Formation in AT1 area，Tahe Oilfield［J］.Lithologic Reservoirs，2010，22（3）：53-58.

［2］付憲弟，王勝男，張亮，等.松遼盆地榆樹林—肇州地區葡萄花油層沉積相類型及沉積演化特征［J］.巖性油氣藏，2013，25（2）：26-30. Fu Xiandi，Wang Shengnan，Zhang Liang，et al.Sedimentary types and evolution characteristics of Putaohua reservoir in Yushulin-Zhaozhou area，Songliao Basin［J］.Lithologic Reservoirs，2013，25（2）：26-30.

［3］王寧，何幼斌，王思文，等.惠民凹陷商河油田商三區沙三上亞段濁積扇沉積特征［J］.巖性油氣藏，2014，26（3）：38-44. Wang Ning，He Youbing，Wang Siwen，et al.Sedimentary characteristics of turbidite fan in upper submember of Es3in the third area of Shanghe Oilfield，Huimin Sag［J］.Lithologic Reservoirs，2014，26（3）：38-44.

［4］GonzalezRC，WoodsRE.Diaitalimageprocessing［M］.2ndEdition. Beijing：Publishing House of Electrics Industry，2002：115-116.

［5］李慶揚，王能超，易大義.數值計算［M］.第5版.北京：清華大學出版社，2008：98-101. Li Qingyang，Wang Nengchao，Yi Dayi.Numerical analysis［M］. 5th Edition.Beijing：Tsinghua University Press，2008：98-101.

［6］宋來忠，王乾峰，彭剛，等.具有大量橢圓隨機分布區域的數值模擬及應用［J］.計算力學學報，2012，29（5）：662-667. Song Laizhong，Wang Qianfeng，Peng Gang，et al.Numerical simulation for the domain with large number of random ellipse grains and their application［J］.Chinese Journal of Computational Mechanics，2012，29（5）：662-667.

［7］Rubin D M.A simple autocorrelation algorithm for determining grain size from digital images of sediment［J］.Journal of Sedimentary Research，2004，74（5）：160-165.

［8］Buscombe D.Estimation of grain size-distributions and associated parameters from digital images of sediment［J］.Sedimentary Geology，2008，210：1-10.

［9］Tedla F K，Tuhtan J A，Hartmann S，et al.Surficial sediment size distribution through the use of digital imaging and autocorrelation analysis［C］.33rdIAHRCongress：WaterEngineeringforaSustainable Environment，2009，496-504.

［10］Warrick J A，Rubin D M，Ruggiero P，et al.Cobble cam：Grain-size measurements of sand to boulder from digital photographs and autocorrelation analyses［J］.Earth Surface Processes and Landforms，2009，34：1811-1821.

（本文編輯：王會玲）

A new method of determining grain size based on rock section image

Yuan Rui1，2，Zhu Rui1，Qu Jianhua3，Sun Yuqiu2，Tang Yong3，Pan Jin1
（1.School of Geosciences，Yangtze University，Wuhan 430100，China；2.School of Information and Mathematics，Yangtze University，Jingzhou 434023，Hubei，China；3.Research Institute of Exploration and Development，PetroChina Xinjiang Oilfield Company，Karamay 834000，Xinjiang，China）

Despite technological advances in lab instruments，grain-size analysis has many limitations，such as low speed. A theoretical method of determining grain size based on spatial autocorrelation coefficient of simulated rock section was proposed.Firstly，spatial autocorrelation coefficient was obtained from a group of known distribution grain size of rock section.Secondly，an unknown distribution grain size of rock section was used to calculate spatial autocorrelation coefficient.Finally，linear least squares method about spatial autocorrelation coefficient was solved with constrains.In order to show the feasibility and availability of this method，a serial theoretical rock sections were simulated by random ellipse process.Relationship between spatial autocorrelation coefficient and grain size of simulated rock section was analyzed.With the decrease of offset or increase of grain size of rock section，spatial autocorrelation coefficient is increasing.Grain size distribution of simulated rock section was determined accurately.For example，gritstone was separated into 0.5～1.0 mm，1.0～1.5 mm and 1.5～2.0 mm，whose computed percentages are respectively 55.8%，24.6% and 20.2%by the proposed method，closing to the actual values，and the variation trend is same as the actual.

spatial autocorrelation coefficient；grain size analysis；least squares method；simulated rock section

P583

A

1673-8926（2015）05-0104-04

2015-03-04；

2015-05-20

國家自然科學基金“珠江口盆地漸新世—中新世潮汐沉積與潮汐周期研究”（編號：41302096）和湖北省自然科學基金“隨鉆跟蹤海量數據處理與成圖關鍵技術研究”（編號：2013CFA053）聯合資助

袁瑞（1987-），男，長江大學在讀博士研究生，研究方向為石油測井與地質。地址：（430100）湖北省武漢市蔡甸區大學路111號長江大學地球科學學院。E-mail：yuanrui87@163.com。