岩石电镜扫描图像的分形特征研究

左 婧 徐卫亚 王环玲 孟庆祥

（1.河海大学 岩土工程科学研究所，南京 210098；2.河海大学 岩土力学与堤坝工程教育部重点实验室，南京 210098）

电镜扫描技术对于分析岩石的微细观结构，进而揭示岩石的力学机理有着重要的作用.随着扫描电镜与计算机图像处理技术的发展，通过考察岩石断面的微观形貌特征，可以在一定程度上反映岩石的力学行为.朱珍德等［1］取锦屏二级水电站引水隧洞大理岩分别进行高水压、高围压、低围压作用下全应力-应变过程三轴压缩对比试验，然后，对大理岩破坏断裂断口进行微观电镜扫描试验，分析不同工况条件下大理岩断口微观形貌特征.郝宪杰等［2］基于电镜扫描实验，对导流洞内柱状节理玄武岩的卸荷破坏机制进行研究.此外，很多学者已经针对岩石流变破坏［3］、岩爆［4］、冲击倾向煤［5］、砂岩断口［6］、大理岩断口［7］等方面进行了相关的研究.因此，利用SEM实验所获得的岩石断口几何图像判定岩石的力学性质，建立岩石微观形态与岩石力学特性的联系，是非常有意义的.

但是目前对于电镜扫描图像的描述多是定性的，岩石作为一种天然材料结构包含了大量不同大小不同层次的孔隙和微裂纹，具有分形结构.分形理论作为近年来在非线性科学中发展出来的新型概念，可用于描述形状复杂、不规则的物体，目前已经在岩土工程中有了一定的应用.谢和平［8］基于前人的成果较为系统地研究了分形在岩石力学中的应用，提出了空隙分形维数的计算方法.Uthayakumar等［9］采用分形理论用来观察酸化后的土壤孔隙变化.毛灵涛等［10］以某高速公路路基加固为例，比较了加固前后土体孔隙微观结构分形维数的变化.但是，现实世界的物体可能具有相同的分形维数但是具有不同的分布特征，因此单一靠分形维数是不够的.景观生态学、医学等领域多引入分形理论的另一个重要概念—间隙度指数（Lacunarity Index）来描述空间分布特征［11］.基于上述考虑，本文采用Matlab开发了岩石电镜扫描图像分形维数与间隙度指数的计算程序，选取了来自某地下洞库风化后花岗片麻岩和煌斑岩脉，通过岩石电镜扫描图片进行计算，探讨了图片放大倍数及阈值对分形维数、间隙度计算的影响，分析了岩石表面结构与分形特征之间的关系与规律.作为一种定量的描述，本文对于相关领域的研究具有一定的意义.

1 图像处理及分形维数、间隙度指数原理与程序设计

1.1 分形维数计算方法

分形维数的测定方法有很多种，如改变观察尺度求维数，根据测度关系求维数，用分布函数求维数，用频谱求维数等.本文采用盒维法来计算电镜扫描图像的分形维数.

盒维数法的基本思路：取边长为r的盒子，将需要测定分形维数的图片覆盖起来，盒子之间不重叠.由于扫描图片已经过二值化，由黑白两部分组成，记录含有黑色部分的盒子数标记为N（r），N（r）的大小取决于盒子边长r.缩小盒子的尺寸r，重复之前的过程，所得到的N（r）就会增加.当r趋向0时，得到盒维法的分形维数：

公式（1）中计算得到的D值即为盒维法测得的分形维数，但在实际应用中只能取有限的边长r.通常根据不同尺寸的r及其对应的N（r），利用最小二乘法拟合以下方程：

公式（2）中计算得到的直线斜率Ds就是盒维法所计算的到的分形维数值.

1.2 间隙度指数计算方法

间隙度指数的概念最早是1983年由Mandelbrot提出的，主要是为了研究间隙的分布［12］.Allain和Cloritre在1991年采用滑动格子法则进行孔隙度指数的计算，并用于分析景观的固定或随机的分维，研究景观的空间格局［13］.滑动格子算法的基本思路是将所研究的区域划分成为网格，记录网格中所研究对象出现或缺失的情况.以不同边长的滑箱从采样区的左上角向右或向下滑动，保证采样点间有重叠部分.记录每次采样时滑箱内研究对象出现的频数，当滑箱滑动过整个区域后，对数据进行统计，计算在不同尺度下的间隙度指数.

间隙度指数的具体计算公式［11］为

公式（3）～（5）中，Λ（r）为孔隙度指数；N（r）为以r×r滑箱滑动整个区域采样的总次数，其计算公式为N（r）＝（M／r＋1）2，M 为区域的边长；S 为景观组分在滑箱中出现的频数；n（s，r）为出现s个景观组分的滑箱个数.

1.3 程序设计

由SEM实验所得到的图片一般为灰度图，无法直接进行分形维数和间隙度指数的计算，在此之前需要对图片进行预处理：圆滑、去噪、二值化，将图像文件转化为0和1的矩阵，然后对所得到的二值化图像进行分形维数和间隙度的计算.根据1.1和1.2节所述的基本原理，本文分形维数、间隙度指数计算流程如图1所示.

图1 分形维数与间隙度指数程序设计

1.4 程序验证

为了验证本文程序的正确性，本文采用Koch曲线作为验证图像，将本文程序的计算结果和景观生态学软件Apack的计算结果进行对比，结果见表1.

虽然Koch曲线分形维数的理论值为1.2618，但是考虑到计算机分辨率有限，无法完全呈现Koch曲线的结构，因此计算结果都小于理论值［14］.对照本程序和Apack软件的计算结果可见采用本程序计算分形参数是正确的.

表1 本文程序与Apack计算结果对比

2 电镜扫描实验

2.1 实验仪器

扫描电镜是一种新型的电子光学仪器，它具有制样简单、所得图像分辨率高等特点，数十年来已被大量应用于生物、医学等学科领域当中，近几年来也被用于观察金属等材料的表面特征.

在实验中，使用了JSM-5610LV型扫描电镜，分辨率在高真空下可达到3.5nm，低真空可达到4.5 nm；放大倍率在18～30万倍之间；电镜配有高灵敏度的二次电子探头和背散射探头，可对固态样品表面形貌进行高速观测，同时配备能谱仪，可在5B～94Pu范围内对样品表面的化学成分组成进行定性以及半定量分析.

图2 JSM-5610LV型扫描电镜

2.2 实验过程及结果

取边长约为1cm某地地下水封洞库风化后的花岗片麻岩和煌斑岩，进行打磨烘干，放入镀膜机内镀膜，之后将其放入载物台并用导电胶将其固定并装入机箱进行照片采集.得到的照片如图3所示.

图3 扫描电镜下分别得到的风化后花岗片麻岩和煌斑岩的表面图

从SEM图片可以发现风化后的花岗片麻岩表面呈鳞状，表面极为不平整，孔隙发育；相比之下煌斑岩的表面平整得多，孔隙不易辨别.

3 结果分析

3.1 图片阈值、放大倍数对分形维数的影响

3.1.1 图片阈值设定对分形维数计算结果的影响

将扫描电镜得到的灰度图转化为二值图时，需要对灰度图进行阈值化.二值化后的图像中，黑色空间代表孔隙结构，白色空间代表骨架结构颗粒.阈值的不同，会导致图片中孔结构和骨架结构的不同，在计算分形维数时的结果就会不同.

选取风化后花岗片麻岩放大倍数250倍下的扫描电镜图片为例，通过Photoshop人为设定阈值，阈值在0～255之间，选取阈值分别为40、50、70、80、100、128、160、200时二值化图片并通过盒维法计算得到相应分形维数结果，如图4（图中同时标出具有代表性的3个不同阈值所得的二值化图片）.

图4 相同风化后花岗片麻岩SEM图片在不同阈值下分形维数计算结果

随着阈值的增大，分形维数先增加后减少，存在最大值.在图片分别为40和140的时候，图片表面形态与真实图片相差过大，无法反映出正常颗粒之间的联通状况，计算所得到的分形维数也不准确.而阈值为70～90区间的时候，较符合真实情况.研究表明当二值化后的图片表面形态与真实图片差别较小时，分形维数都比较大.由此可以发现，图像质量随着阈值的变化而显著变化，同时对分形维数的计算影响巨大，因此最佳阈值的选取对分形维数具有重要影响.3.1.2 放大倍数对分形维数计算的影响

选取花岗片麻岩在放大倍数分别为250、500、1000、2000、5000、10000情况下的电镜扫描照片（放大倍数大的图片区域从放大倍数小的图片区域中选取），为避免阈值的干扰，将二值化阈值设置为相同值，使用盒维法计算分形维数.

增加图片放大倍数时，扫描范围更小，所观察到的孔隙结构更加详细，图片表面结构变化复杂程度变低，骨架结构对孔隙填充能力降低，因此导致分形维数结果随着放大倍数增加而减小.结果如图5所示.

图5 花岗片麻岩图片不同放大倍数下分形维数变化曲线

根据图5所示的分形维数变化曲线，可见随着放大倍数的增加，当放大倍数大于2000后，分形维数基本稳定.因此，在研究图像的分形特征上应该选用分辨率较高的电镜扫面图像.

3.2 图像阈值、放大倍数对间隙度指数的影响

3.2.1 图片阈值设定对间隙度计算结果的影响

正如上文所提到的，阈值的不同改变图像中孔结构和骨架结构的差异，阈值设置越高，黑白图中黑色部分越多，相对白色成分越少；即孔结构越多，骨架结构越少.同样选取风化后花岗片麻岩250放大倍数下的扫描电镜图片为例，设定不同阈值，运用滑箱运算法则进行间隙度计算.

图6 风化后花岗片麻岩SEM图片不同阈值下的孔隙间隙度计算结果图

进行孔隙间隙度指数计算时，固定滑箱变长，阈值设置越高，代表孔隙结构的黑色区域覆盖率越高，出现频率越高，相应间隙度指数越小，如图6（a）所示.选取滑箱边长为横坐标，可以看到图6（b）中不同阈值下滑箱边长与间隙度指数曲线的变化情况.阈值越小，曲线越远离横坐标轴，阈值越大，反之.通过间隙度指数在景观生态学中的应用可以了解到3种景观分布特征［11］.阈值越大，孔隙间隙度指数越小，表明孔隙分布越均匀，但同时也无法体现电镜扫描图像的分布特征.

综上所述，计算间隙度指数时阈值的选取至关重要，只有正确的阈值才能反映出合理的孔隙结构、骨架结构和研究对象的分布性，得到较为准确的间隙度指数.

3.2.2 放大倍数对间隙度指数及孔隙分布的影响

放大倍数对间隙度指数的影响如图7所示，可见随着放大倍数的不同，岩石电镜扫描图像的孔隙分布可以分为两种类型.当放大倍数为250和500时，由于扫描图片的范围较大，孔隙分布表现出均匀分布的特征.当放大倍数为1000和2000时，间隙度指数曲线基本相同.这是由于随着放大倍数的增加，扫描图片选取范围越来越小，间隙度指数越高，曲线离横坐标越来越远，孔隙分布越不均匀，孔隙分布体现异质性，能够较为准确地反映岩石的细观结构.

由此可见，放大倍数的大小对孔隙分布是否均匀的判断有所影响.因此在计算间隙度指数的时候尽量选取放大倍数大的图片，这样可以大范围内地展现出研究对象的结构分布，真实反映出分布情况.

图7 放大倍数对间隙度指数的影响

3.3 两种岩石SEM图片分形维数计算、间隙度指数对比分析

3.3.1 两种岩石SEM图片分形维数计算对比分析

同时选取了相同放大倍数的电镜扫描图片采用盒维法进行分形维数的计算统计，阈值采用最佳值，计算结果见表2.

表2 各个放大倍数岩样的分形维数值

通过表2可以发现，相同放大倍数情况下，风化后的花岗片麻岩的分形维数比煌斑岩脉的大，可以推断岩石表面风化程度和孔隙发育程度可以用分形维数定量表明，岩石表面越粗糙，风化和孔隙发育程度越高，分形维数相对越高.

3.3.2 两种岩石SEM图片间隙度指数对比分析

选取最佳阈值、放大倍数2000下的图片进行两种岩石孔隙间隙度指数的分析.相同滑箱边长时，风化后花岗片麻岩间隙度指数比煌斑岩脉的值小，而随着滑箱边长减小，风化后花岗片麻岩孔隙间隙度指数变化程度比煌斑岩脉的程度小，表明风化后花岗片麻岩表面空间格局较煌斑岩脉更为复杂.结果见图8.

图8 两种岩石孔隙间隙度指数对比

4 结 论

通过岩石电镜扫描实验及其分形特征分析，主要结论如下：

1）电镜扫描图像在进行二值化时，随着阈值的增大，盒维法计算所得到的分形维数先增加后减少；孔隙的间隙度越小，孔隙分布越稠密，异质性特征差异越不明显.因此在计算分形维数时应该选用分形维数最大值对应的阈值作为最优阈值，使得二值化后的图片更加真实反映了岩石颗粒表面的形态.

2）设定相同阈值时，随着放大倍数增加，分形维数结果逐步降低，孔隙间隙度指数越大，孔隙分布表现越不均匀，由此可见对于低渗透岩石，电镜扫描图像放大倍数至少大于1000倍才能体现出岩石的分形结构特征.

3）选取最佳阈值和放大倍数风化后的花岗片麻岩的分形维数比煌斑岩脉的大，岩石表面越粗糙，风化程度和孔隙发育程度越高，分形维数相对越高.随着电镜扫描在孔隙结构研究方面的研究和应用越来越广泛，分形理论作为一种定量描述方法，对于揭示岩石的力学机理有着重要的意义.

［1］朱珍德，张 勇，徐卫亚，等.高围压高水压条件下大理岩断口微观机理分析与试验研究［J］.岩石力学与工程学报，2005，24（1）：44-51.

［2］郝宪杰，冯夏庭，江 权，等.基于电镜扫描实验的柱状节理隧洞卸荷破坏机制研究［J］.岩石力学与工程学报，2013，32（8）：1647-1654.

［3］张强勇，杨文东，陈 芳，等.硬脆性岩石的流变长期强度及细观破裂机制分析研究［J］.岩土工程学报，2011，33（12）：1910-1918.

［4］谭以安.岩爆岩石断口扫描电镜分析及岩爆渐进破坏过程［J］.电子显微学报，1989，8（2）：41-48.

［5］宋晓艳，王恩元，刘贞堂，等.冲击倾向煤破坏断口微观特征及其机理研究［J］.中国矿业大学学报，2008，37（6）：775-779.

［6］左建平，谢和平，周宏伟，等.温度-拉应力共同作用下砂岩破坏的断口形貌［J］.岩石力学与工程学报，2007，26（12）：2444-2457.

［7］黄润秋，黄 达.高地应力条件下卸荷速率对锦屏大理岩力学特性影响规律试验研究［J］.岩石力学与工程学报，2010，29（1）：21-33.

［8］谢和平.分形几何及其在岩土力学中的应用［J］.岩土工程学报，1992，14（1）：14-24.

［9］Uthayakumar R，Prabakar G A，Azis S A.Fractal Analysis of Soil Pore Variability with Two Dimensional Binary Images［J］.Fractals，2011，19（4）：401-406.

［10］毛灵涛，薛 茹，袁则循.软土路基微结构扫描电镜图像的分形［J］.长安大学学报：自然科学版，2007，27（2）：30-33.

［11］Plotnick R E，Gardner R H，O'Neill R V.Lacunarity Indices as Measures of Landscape Texture［J］.Landscape Ecology，1993，8（3）：201-211.

［12］Bunde A，Havlin S.Fractals and Disordered Systems［M］.Springer-Verlag New York，Inc.1991.

［13］Allain C，Cloitre M.Characterizing the Lacunarity of Random and Deterministic Fractal Sets［J］.Physical review A，1991，44（6）：3552.

［14］李水根.分形［M］.北京：高等教育出版社，2004.