张 欣，张 栋，杨 倬，高世臣，张 艳

(1. 中国地质大学(北京) 数理学院，北京 100083； 2. 中国石油长庆油田分公司 勘探开发研究院，陕西 西安 710018；3. 中国地质大学(北京) 地球物理与信息技术学院，北京 100083)

利用岩石薄片进行镜下岩石颗粒特征分析、岩石成分鉴定、孔隙类型识别及其发育程度的表征是复杂岩性油气藏沉积特征、成岩作用以及储层综合评价等研究工作中极其重要的技术手段(Moorhouse，1986； 林培英，2005)。近年来随着岩性油气藏勘探开发工作的深入发展，大量的取心资料一方面需要快速地完成鉴定分析工作，另一方面沉积、成岩极为复杂的致密储层评价对岩石薄片鉴定在定量化、精确化方面又提出了更高的要求(张厚福等，2008； 李富恒等，2009)。传统的岩石薄片鉴定通常是专业人员交互式的手工操作进行识别和统计，以定性-半定量为主，且耗时较多。自20世纪60年代至今，数字图像处理领域一直在发展与推进，已在医学和空间项目的基础上用于更广泛的范围，并且在地质领域的岩相学问题、油气勘探与预测、油气储量计算等方面也发挥了重要作用，早在1993年的油气开发领域中便已得到应用。利用岩石切面或者岩石薄片的图像分析，进而获取岩石参数或进行岩石特征推断也进一步促进了数字图像处理技术在地质学领域的应用。数字图像处理技术的引入使得岩石矿物的镜下鉴定和分析在一定程度上得以简化，工作效率也大为提高(吴拥等，2013； Chauhanetal., 2016)。但是岩性油气藏致密储层的成岩极为复杂，通常压实作用很强，胶结作用严重，原生岩石颗粒之间接触紧密、边缘模糊，而后期溶蚀作用的改造使得部分岩石颗粒轮廓模糊，与孔隙之间的边缘变得复杂。岩石颗粒溶蚀残余、孔隙内的自生填隙物等在镜下薄片图像中形成较多的噪声点，而传统的边缘提取算子、阈值分割技术大多利用灰度信息及邻域信息进行边缘提取，在面对高信息量、高复杂度的数据时具有局限性，岩石颗粒分割和孔隙提取效果较差。本文提出了一种新的岩石薄片图像采集和数字化处理分析方法，以期改善岩性油气藏复杂储层岩石薄片颗粒分割和孔隙提取效果，提高复杂岩石薄片的鉴定效率和可靠性。

1 角度域正交偏光序列图采集方法

岩石薄片的制作通常是将垂直岩石层理方向的小块岩样磨制成厚度约0.03 mm的薄片，固定于载玻片上，并覆盖盖玻片。通过光学显微镜观察可鉴定分析岩石薄片中的矿物、孔隙等特征(叶润青等，2011； Freezetal., 2012)。在单偏光下观察时，光线通过位于载物台下的下偏光镜，射入薄片样品的平面偏光为单一振动方向，多数透明矿物颗粒无明显的颜色显示，岩石颗粒之间的边缘特征不显著；如果增加位于物镜之上的上偏光镜，在正交偏光观测模式下，岩石矿物颗粒显示出明显的光学特征差异，例如矿物颗粒的干涉色、纹理特征、颗粒边缘凸起差异等。从晶体光学原理分析，这种差异性取决于矿物的成分、颗粒方位等因素，通过旋转载物台，连续改变岩石矿物颗粒光率体方位与正交光之间的夹角，可以跟踪岩石颗粒的消光变化特征，进而综合多角度位置上的特征判别颗粒边缘及其矿物成分，以提高岩石薄片的鉴定可靠性。通过与显微镜连接的照相机进行多角度正交偏光的岩石薄片图像采集，进而采用图像数字化处理技术进行岩石颗粒分割、粒度分析、矿物成分鉴定和孔隙提取等储层特征综合研究工作。

测试岩石样品为鄂尔多斯盆地某区长6段砂岩薄片，前期人工鉴定主要岩石颗粒成分有石英(26.43%)、长石(44.31%)及岩屑(29.26%)，面孔率24.20%。采取连续偏光图像采集的方法，以10°为间隔连续旋转载物台，进行360°图像采集，共得到36张正交偏光序列图像，根据矿物在正交光镜下的消光机理可知，消光周期为90°(Moorhouse，1986； 林培英，2005)，因此36张图像为4个完整的消光周期(图1)。图像采集过程中，显微镜放大倍数为4倍，通过图像角度方位校正对齐处理，获取的图像像素为1 000×1 000的圆形视域。

2 岩石矿物颗粒角度域消光特征分析

选取典型的石英、长石、火山碎屑、孔隙各2块，对目标颗粒0°～360°偏光角度范围内的同一视域进行标定，视域中心点位置坐标对应图1中20°偏光图中的标注点。计算各偏光角度的视域内灰度均值，形成矿物的角度域灰度特征曲线，统一曲线的初始相位，绘制3种类型矿物灰度的角度域特征图谱(如图2a所示)，并进一步对矿物颗粒图像R、G、B三通道特征在角度域的变化特征进行分析(图2b～2e)。

图 1 正交偏光显微镜下0°～90°周期内岩石薄片序列图Fig. 1 Multi-angle (0°～90°) microscopic images of thin section (crossed nicols)Q—石英； F—长石； V—火山碎屑； P—孔隙Q—quartz； F—feldspar； V—volcanic debris； P—pore

颗粒灰度存在角度域的周期性变化，且以90°为一个消光周期，周期的角度域信息完整表达了不同矿物颗粒的特征。对目标薄片的3种矿物颗粒多次选择并统计，普遍情况下，长石颗粒的灰度值低于石英颗粒，即镜下光度较大，岩石颗粒灰度值相对大小顺序为：石英>长石>火山碎屑，石英曲线的峰值主要分布在200～255之间，长石曲线的峰值主要分布在75～150之间，火山碎屑的灰度最低，其峰值主要分布在25～65之间。除不同类型颗粒间灰度差异外，颗粒R、G、B通道值存在明显不同的相互关系。石英颗粒的R通道峰值相比于G、B两通道明显较低，G、B通道峰值接近但B通道振幅略小；长石与火山碎屑的三通道值差距较大，由高到低按照B、G、R依次排列，其中长石颗粒的G、R通道差异大，而火山碎屑的B、G通道差异大。岩石颗粒R、G、B三通道特征显示了不同矿物颗粒色调的差异，而对于同一颗粒而言，随着消光角度的变化，其光度呈现有规律的明暗变化，但其色调始终是保持不变的。岩石颗粒正交偏光镜下特征在角度域的变化、不同岩石颗粒变化特征的差异性是进行岩石薄片数字化鉴定分析的基础。

3 正交偏光图像岩石颗粒分割和孔隙提取

3.1 常规方法提取岩石颗粒边缘

岩石薄片序列图中岩石颗粒表面纹理复杂，颗粒内部以及边缘附近的噪点较多，严重影响了边缘提取的效果，造成边缘零散，提取结果较粗，岩石颗粒形态、孔隙结构和岩石颗粒溶蚀特征等难以显现。双边滤波在消除颗粒内部及边缘噪点的同时，对颗粒边缘特征及孔隙特征弱化程度较低，多数颗粒边缘更加清晰，因此选择双边滤波的去噪方法对岩石薄片序列薄片进行了预处理。

图 2 不同岩石组分角度域正交偏光消光谱特征Fig. 2 Spectroscopies in angular domains of different components under crossed nicols microscopya—各岩石组分灰度； b—石英颗粒； c—长石颗粒； d—火山碎屑； e—孔隙内部a—gray level of component； b—quartz； c—feldspar； d—volcanic debris； e—pore

3.1.1 Sobel、Canny算子的单偏光角度边缘检测

采用常用的图像边缘提取方法Sobel算子和Canny算子进行薄片岩石颗粒边缘进行提取。图3是薄片在20°和70°角度时正交光镜下的岩石颗粒边缘及轮廓特征图。

传统的岩石薄片边缘提取技术，多基于岩石薄片在单一角度位置的正交光镜下图像，提取结果通常无法准确提取全部岩石颗粒的轮廓。从图3可以看出，几乎所有矿物颗粒的轮廓特征都会在某些消光角度的图像中有所显现，但是同一颗粒的轮廓特征极少在所有消光角度的图像中均有明显显现。因此，采用单一角度正交偏光图像对岩石薄片进行轮廓提取存在一定的缺陷：单一角度正交偏光图像中总会有部分岩石颗粒的边缘特征不明显，难以提取完整的边缘信息，颗粒划分粗糙。对于某一角度正交偏光镜下，两相邻颗粒光度接近时，常规方法很难精确提取颗粒边缘和轮廓。

3.1.2 基于角度域序列图像的边缘融合

针对单一角度正交偏光图像只能较好地刻画出部分岩石颗粒的边缘特征，且某一岩石颗粒的边缘总会在某些角度的图像中有明显显现的特点，采用角度域序列图像边缘提取结果叠加的方式进行边缘特征的融合。图4为岩石薄片在0°～90°一个消光周期内的10张图像R、G、B三通道参数边缘提取结果各自的叠加图以及三通道叠加融合图。

图 3 基于灰度的单偏光图像边缘提取Fig. 3 The edge diagram extracted in gray level of plainlight imagea—Sobel，20°； b—Sobel，70°； c—Canny，20°； d—Canny，70°a—Sobel at 20°； b—Sobel at 70°； c—Canny at 20°； d—Canny at 70°

相比于传统的单一角度的正交光图像边缘提取的结果，基于角度域和颜色通道的边缘提取和叠加融合方式实现了颜色通道信息和消光角度的纵横信息联合拓展。对比Sobel和Canny两种边缘提取方法在R、G、B和灰度参数中提取边缘的叠加融合结果，可以看出不同通道参数展现的颗粒边缘特征较为一致，而多角度序列图像的叠加，展现了全部颗粒的边缘和轮廓，且较为清晰。但是这种多图像边缘的直接叠加也导致多图像噪点的叠加，给岩石颗粒的准确分割带来困难。

图 4 Sobel及Canny算子周期内三通道的边缘提取与叠加Fig. 4 The overlay edge maps of three channels in a cycle with Sobel and Canny operatorsa～d—Soble算子边缘； e～h—Canny算子边缘； 按列依次为R、G、B通道，最后一列为三通道叠加结果a～d—the edge diagrams with Sobel operator； e～h—the edge diagrams with Canny operator； each column is followed by R, G, B, and three channel overlay results

3.2 基于角度域消光特征的相关分析方法

3.2.1 序列消光图像相关分析方法

简单地对不同角度的边缘提取结果进行叠加，是对各角度边缘提取结果取并集的结果，当各角度边缘图参数设置不佳时，对最终提取结果有较大影响，若保留较多单角度边缘信息，会造成最终叠加图的过分割。该方法虽然在最大程度将边缘复刻至最终结果图，但并未充分考虑各角度域图像之间的内在联系，无法表达岩石颗粒角度域消光特征的变化规律，无疑是一种有效信息的损失。因此，构建像素点一个消光周期内的灰度向量A(x,y)作为该像素点的完整消光特征表达：

A(x,y)=vec(f0(x,y),f1(x,y),…f9(x,y))

其中(x,y)表示像素点位置，fi(x,y)表示(x,y)在第i幅序列图的灰度值。

以Pearson相关系数来度量像素点之间的相似性，选定目标像素点，计算目标像素点与其邻域内的各像素点灰度向量间的相关系数。图5是对图像不同颗粒和边缘101×101像素邻域相关系数的计算结果，图5a～5f分别对应图1中70°偏光图像中的标注区域。图中越接近红色的区域与中心目标点融合效果越好。颜色相近的区域为强相似区域，可作为融合后的分区，一般表征同一颗粒、边缘或噪声点的内部。不同颗粒和不同边缘组合位置的相关性特征如下：

(1) 同一颗粒内部的点之间相关系数大(大于0.75)，呈现明显的独立分区，且岩石颗粒分割特征显著(图5a)；

(2) 边缘区域与颗粒内部在相关系数上存在明显差异，相邻颗粒形成明显的分区，而颗粒之间边缘附近的相关系数则明显较低(小于0.3)；低亮度的岩石颗粒边缘特征偏弱，分区特征不明显(图5b、5c)；

(3) 岩石颗粒之间的孔隙内部相关特征明显，孔隙和颗粒之间存在渐变特点(图5d、5e)；

(4) 溶蚀的长石颗粒边缘轮廓不明显，内部结构复杂，溶蚀颗粒残余部分和完全溶蚀成为孔隙的部分相关系数较高，而孔隙填物分布区相关系数较低甚至负相关(图5f)。

3.2.2 基于角度域相关的边缘提取方法

分别选取石英颗粒内部、长石颗粒内部、孔隙结构内部及边缘点作为目标点，图6为边缘提取示意图，图6b中的每一个黑色方框区域即为指定的窗口大小。以采用5×5窗口大小的图6a为例，其中的数字5为目标点，9个数字分别对应于图7中每个图形的9条线，最终产生图7a～7d中4种情形下9点的消光特征雷达图。

边缘点作为相邻两颗粒的公共边，所选邻域跨越两颗粒或一颗粒，包含颗粒内部点，因此部分点与目标点相似程度差，为表征这种特点，采用目标点与邻域内其他点的相关系数均值、标准差作为边缘提取的评判指标，并进一步构建相关系数均差比指标，方法如下：

图 5 不同颗粒和边缘组合位置的邻域相关性Fig. 5 Neighborhood correlation of different grains and edge patternsa—颗粒内部； b—相交边缘； c—单侧低亮边缘； d、e—粒间孔隙； f—长石溶蚀区域a—inside grain； b—intersection edge； c—gray edge on one side； d, e—intergranular pores； f—feldspar dissolution zone

图 6 基于相关的边缘提取示意图Fig. 6 Diagram of edge extraction based on correction

向量之间的相关系数计算公式为：

r(A(x,y),A(m,n))=Cov(A(x,y),A(m,n))/

σ(A(x,y))·σ(A(m,n))

目标点相关系数均值：

目标点相关系数标准差：

目标点相关系数均差比指标：

rms(x,y)=rmean(x,y)/rstd(x,y)

其中(x,y)为目标点，δ为正方形窗口，k为窗口大小，(m,n)为窗口内非目标点。

由图7可知，石英和长石内部均表现为强相关；孔隙内部相关性有扰动，为非均质结构；边缘相关性变化最大，表现为相关均值较低，相关标准差较高，均差比明显小于颗粒内部。

图 7 5×5邻域9点的消光特征雷达图Fig. 7 Radiation diagrams of extinction spectroscopies at 9 points in 5×5 neighborhoodsa—石英内部； b—长石内部； c—孔隙内部； d—边缘a—quartz； b—feldspar； c—pore； d—edge

3.2.3 基于角度域相关的边缘提取及岩石颗粒分割和孔隙提取

采用上述方法，对岩石薄片在0°～90°一个消光周期内的10张图像进行相关性分析，并提取相关系数均值和相关系数标准差(图8a、8b)。对比于传统的Canny、Sobel算法，基于相关分析的边缘提取结果轮廓的闭合性更好，弱边缘的表达更为清晰，孔隙结构更为明显。

利用相关系数均值图对岩石颗粒分割和孔隙进行提取。其中岩石颗粒和孔隙内部的相关系数均大于0.75，孔隙部分在正交偏光镜下图像的灰度值低于25，且R、G、B三通道参数相近，颗粒与孔隙的边缘的相关系数均值小于0.3。

图 8 岩石薄片颗粒分割成果图Fig. 8 Grain segmentation results of thin sectiona—相关系数均值； b—相关系数标准差； c—岩石颗粒分割图； d—孔隙分布图a—mean coefficient of correlation； b—correlation coefficient standard deviation； c—grain segmentation map； d—pore distribution map

3.2.4 效果评价

从边缘提取效果层面分析，比照图3，相比于传统的单偏光方法，从多偏光序列图像中提取的边缘信息更加完整；相比于图4的Canny、Sobel算子，该方法对边缘的粗细刻画更为准确，精确到像素级，且有效降低了图4中算法本身造成的毛细边缘太多而产生的过分割现象，为进一步的岩石矿物成分识别、定量分析及参数统计打下了良好的基础。

从定量分析效果层面分析，统计图8c和图8d岩石颗粒和孔隙分割图像中孔隙部分所占比例，即面孔率，为25.36%；采用人工镜下网格交点法统计的面孔率为24.20%，两者偏差为1.16%，总体结果较为一致。网格法作为基于人为经验判断边缘进而估计所占区域大小的传统方法，其统计精度约为网格大小的一半，对微小孔隙、溶蚀颗粒模糊边缘等区域的统计易产生较大误差，在精细定量化方面传统方法略显薄弱。另一方面，多偏光序列图像相关分析法尽管在统计精度、定量化和自动化等方面优点突出，但是对岩石成分和孔隙类型的准确人工鉴定依然是对该方法分析结果进行标定和校正的基础。通过两种方法的有机结合和相互印证，可以获得更为准确、更加精确的结果。

4 结论与建议

(1) 岩石薄片角度域序列正交偏光镜下图像的数字分析技术，一方面可以提取并分析岩石矿物在完整的消光周期内的连续变化特征，更为全面地揭示出岩石矿物的光学特征，另一方面岩石薄片图像的RGB三通道分析则可以提供岩石矿物在颜色通道维度上的变化特征，很好地改善了岩石颗粒边缘特征的刻画效果、矿物成分鉴定和识别效果，并提高了工作效率。

(2) 综合角度域正交偏光消光特征的岩石薄片图像相关分析法，综合了岩石薄片图像的角度域消光变化特征和像素点邻域相关性，与Soble和Canny等边缘提取算法相比，明显提升了岩石颗粒边缘的提取效果，更为准确地刻画了岩石颗粒的轮廓和孔隙分布特征。

(3) 岩石薄片角度域正交偏光消光谱相关性分析表明，在同一矿物颗粒和孔隙内部的消光谱相似性较强，邻域相关性较高，多大于0.75，而在矿物颗粒之间的接触边缘、颗粒与孔隙相邻部位则相关系数较低，多小于0.3，甚至负相关。以此为据可以较为准确地分割出矿物颗粒和孔隙。

(4) 岩石薄片正交偏光镜下的消光特征在角度域上的连续变化，为研究岩石矿物提供了一个新的维度和视角，岩石矿物在角度域上的消光谱特征不仅仅可以用于刻画岩石颗粒形态、孔隙特征等，同时还提供更为丰富的岩石颗粒形态、岩石矿物成分、孔隙类型及其内在填隙物等信息(两个岩石颗粒形态)，为岩石薄片的自动化、数字化、定量化的鉴定分析提供了数据和信息基础。岩石薄片正交偏光消光谱的研究及应用，在沉积岩石粒度和产状分析、矿物成分定量识别、孔隙结构参数提取、孔隙填隙物识别及岩石沉积、成岩特征研究等方面必将发挥极为重要的作用。

(5) 在对岩石薄片的边缘提取的研究中，颗粒边缘的各向异性是一个显著的信息，且边缘和颗粒的敏感窗口大小相距较大，即关于边缘的方向信息与窗口大小对计算效率的提升在本文的研究角度和技术方法中并未探讨，基于岩石薄片信息利用的完整性的探讨和研究仍需进一步展开。