一种基于岩石薄片图像的粒度分析新方法

袁　瑞，朱　锐，瞿建华，孙玉秋，唐　勇，潘　进

（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。