摘要：基于分形理论，利用MATLAB开发密封环端面激光共聚焦扫描图像分形维数计算及分析程序，利用Sierpinski地毯对程序的计算结果进行标定，计算不同分割阈值下的图像分形维数，给出一种分割阈值的选取方法。研究结果表明：设计开发的分形维数计算程序具有较高的准确性；图像分割阈值的选取对密封环端面分形维数的计算有重要影响，给出的分割阈值选取方法，能使图像二值化效果最佳，使扫描图像能最大程度地反映密封环端面形貌的分形特征，从而提高密封环端面形貌定量描述的精确度。
关键词：机械密封；表面形貌；计盒维数；分割阈值
密封环端面的形貌是机械密封加工和工作过程中表现出来的重要特征，它对机械密封的摩擦特性和密封特性具有重要影响[1-2]。研究密封环摩擦副粗糙表面的固有特性，正确表征摩擦副的表面形貌，对了解密封端面间的泄漏通道状况，以及提高机械密封动静环接触端面的磨合质量具有重要的意义[3]。
20世纪70年代，法国数学家Mandelbrot[4-5]正式提出并创立了一种探索自然界复杂形态的数学分支——分形几何学。分形几何能够在更深层次上深刻地描述、研究和分析普遍存在的杂乱无章的、不规则的、随机的自然现象，在摩擦磨损领域中引入分形理论实现了对摩擦副粗糙表面形貌的科学的定量描述，从而推动了机械密封端面形貌与探寻密封泄漏通道的研究[6]。在一些基于分形理论的机械密封端面研究中[7-9]，大多采用W-M或M-B函数来表征和模拟密封端面，同时利用功率谱法和结构函数法来研究轮廓曲线的分形特性，但是函数模拟只能得到统计意义上的形貌，并不能全部反映密封端面的泄漏通道情况，因此利用扫描图像来研究密封端面的形貌特征就具有重要的意义。在以往一些基于扫描图像的分形特征研究中[10-12]，多结合图像处理技术对图像进行预处理，然后利用计盒维数等方式求取二值化后图像的分形维数，然而二值化过程中阈值的选取对图像的特征表现和分形维数计算具有重要的影响[13-14]，所选取的分割阈值应能最大程度地将密封端面扫描图像的分形特征表现出来。
本文作者拟利用密封环端面的激光共聚焦扫描图像，通过MATLAB编程计算密封端面的分形维数，讨论二值化过程中分割阈值选取对分形维数的影响，给出一种图像分割阈值的选取方法，突出图像的自相似特征，从而为密封环端面的分形表征提供一种快速、精确的新方法。
1 基于MATLAB的分形维数算法设计
分形维数是分形几何理论中最重要的基本概念之一，目前已经有很多维数的定义，主要包括Hausdorff维数、计盒维数(CBD)、修正计盒维数、填充维数等。Hausdorff维数是分形几何理论的基础，可以说分形几何的理论体系是建立在这一基础之上，但是Hausdorff维数只适合分形几何的理论推导，只能通过分析的方法获得一小类规则的纯数学分形的Hausdorff维数，它对实际应用中提出的分形维数的计算问题无能为力。鉴于此，人们提出了计盒维数的概念，由于计盒维数物理含义直观，易于进行程序化计算，因此得到了广泛的应用[15]。
1.1 计盒维数的原理
设A是Rn空间的任意非空有界子集，对于任意一个r>0，Nr(A)表示用来覆盖A所需边长为r的n维盒子的最小数目。如果存在一个数DB，使得当r→0时，有
Nr(A)∝1/rDB
(1)
那么称DB为A的计盒维数(简称盒维数)。
进一步求得
(2)
计盒维数这个形式的定义在实际中有广泛的应用。在计算一个平面集A的计盒维数时，可以构造一系列边长为r的正方形，然后计算不同r值的盒子与A相交的个数Nr(A)， 这个维数是当r→0时， Nr(A)增加的对数速率，或者可以由函数lnNr(A)相对于lnr图的斜率的负值来估计[16]。
1.2 程序设计
图像在计算机中是以矩阵的形式存储，一幅长度为M个像素，宽度为N个像素的图像，可以看成是一个M×N的矩阵，矩阵的每一个元素代表一个像素，元素的值是像素点颜色或索引色[17]。利用激光共聚焦获得的密封环端面图像是RGB格式，无法直接进行计盒维数的计算，因此需要对图像进行预处理，包括灰度变换、平滑滤波、降噪等，再选取一个合适的分割阈值，将图像二值化，这样就将图像处理成只包含0(孔隙)和1(固体骨架)的二值矩阵。根据计盒维数的原理，可以通过MATLAB编写程序计算计盒维数，步骤如下：
(1)读取二值化的密封环端面图像，并将其数值矩阵赋值给变量H；
(2)用不同边长尺寸的正方形盒子对矩阵H进行覆盖，统计在不同尺寸r下包含目标像素(0—孔隙)的盒子数，计入数组Nr中；
(3)在双对数坐标中利用最小二乘法对(lnr,lnNr)进行拟合，得到直线斜率的负数即为图像的计盒维数。
程序的流程图如图1所示。
图1 MATLAB计算计盒维数程序流程图
Fig 1 Flowchart of program calculating counting-box dimension using MATLAB
1.3 程序验证
为了验证设计程序的正确性，采用5级Sierpinski地毯作为标定图像，将计算结果与标准结果进行对比，如表1所示。可以看出，本文的计算结果与理论维数相对误差非常小，说明本文作者设计的程序是可行的，同时也具有很高的精确度。
表1 计算的Sierpinski地毯计盒维数与理论维数比较
Table 1 Comparison between the theoretical fractal dimension and the calculated counting-box dimension of Sierpinski carpet
2 实验及分析
2.1 图像获取
采用型号为奥林巴斯OLS4100的激光共聚焦显微镜对材料为硬质合金YG-8的动环以及材料为浸渍酚醛树脂碳石墨的静环端面进行扫描，得到的图像如图2所示。
图2 动、静环激光共聚焦扫描图像
Fig 2 The LSCM images of rotating ring
(a) and stationary ring (b)
可以看出，动环表面孔隙较少，表面较平整，孔隙与骨架易区分；相比之下，静环表面较复杂，孔隙较多，且不易辨别。
2.2 阈值的影响及选取
将激光共聚焦获取的扫描图片进行二值化处理时，需要对灰度图像进行阈值分割，分割后的图像中白色像素代表骨架，黑色像素代表孔隙。阈值的不同，会导致图像中白色和黑色像素的比例、分布结构不同，在计算计盒维数时的结果也会不同。对于一幅8位的灰度图像，每一个像素的灰度范围为0～255，其中0代表黑色，255代表白色，阈值的取值范围为0～255，对应在区间0～1内，按比例选取不同阈值将图像二值化，采样间距为0.05，并计算出图像的分形维数，作出分割阈值与分形维数的关系图，如图3所示。可以看出，随着阈值的增加，分形维数也不断地增加，因此分割阈值的选取，对图像计盒维数的计算结果具有重要影响。
图3 动环扫描图像分割阈值与分形维数的关系
Fig 3 The relation between the segmentation threshold and the fractal dimension of LSCM image for rotating ring
在图像处理领域，通常采用Otsu法(又称为最大类间方差法)来选取图像的分割阈值，它是按图像的灰度特性，将图像分成背景和目标两部分。背景和目标之间的类间方差越大，说明构成图像的两部分的差别越大，当部分目标错分为背景或部分背景错分为目标都会导致两部分差别变小。因此，使类间方差最大的分割意味着错分概率最小。 但是Otsu法对噪声和目标大小十分敏感，仅适合图像目标和背景分布呈正态分布且各像素量和方差基本相当的情况，当目标与背景的大小比例悬殊时，类间方差函数呈现双峰或多峰，此时阈值分割后的图像质量就不理想[18]。这种情况下，就需要通过尝试的办法来选取最佳的分割阈值。
本文作者将阈值的灰度范围0～255等比例对应到0～1区间，并按照0.01的间距进行采样，统计对应比例阈值下孔隙面积的大小，也即图像中黑色像素的个数。如图4所示为动环扫描图像不同分割阈值与孔隙面积的曲线关系图。可以看出，随着阈值的增大，孔隙面积也是不断增大的。对起始段和上升段的数据点进行直线拟合，交点处的阈值约为0.4，同时，当阈值大于0.5时，孔隙面积基本保持稳定，也就是说图像中黑色像素数量达到最大；当阈值小于0.2时，孔隙数目非常少，即图像中基本没有黑色像素存在，在0.2～0.4区间内，曲线变化缓慢，在0.4～0.5区间内，曲线上升幅度明显增加，这说明在拐点之前的图像分析主要集中在孔隙区域，之后则集中在背景区域，因此阈值应当选择图中拐点处对应的阈值作为二值化的最佳阈值[19]。
为了进一步验证此交点处的阈值为最佳阈值，利用MATLAB做出比例阈值为0.2、0.4、0.5时的二值化图像，如图5所示。可以看出，当阈值为0.4时，二值化的图像能较真实地反应孔隙分布，而当阈值为0.2和0.5时，二值化的图像都有严重的失真，不能反映真实的孔隙分布，因此可以选取0.4作为计算分形维数时的最佳阈值。
图4 动环扫描图像分割阈值与孔隙面积关系
Fig 4 The relation between the segmentation threshold and the pore areas
图5 原始图像和不同分割阈值下的二值化图像
Fig 5 The original image (a) and the binarization images of the proportional segmentation threshold of 0.2(b),0.4(c)and 0.5(d)
2.3 分形维数计算
利用最佳阈值对灰度图像进行二值化处理，根据计盒维数的原理，通过MATLAB编写的程序计算分形维数，采用边长为r的正方形盒子对图像的数值矩阵进行覆盖，统计不同边长下，包含目标像素的盒子数，获得的数据如表2所示。
表2 盒子边长r及相应边长下所需的盒子数Nr
Table 2 The side length (r) and the corresponding number (Nr) of boxes
将表中的两组数据分别取对数，并利用最小二乘法在双对数坐标中对其进行拟合，得到的直线的斜率的负数就是该扫描图像的计盒维数。经过计算得到的计盒维数DB=1.797 7，两组数据的相关系数R=-0.997 5，因此可以认为二者具有很高的线性相关性，即密封端面具有高度的自相似特性。计算的结果如图6所示。
图6 盒子边长r及盒子数Nr双对数坐标下的拟合方程
Fig 6 The fitting equation in double logarithmic coordinates of side length (r) and the number (Nr) of boxes
3 结论
(1)利用MATLAB图像处理，结合分形的基本理论，设计开发了基于MATLAB的分形维数计算分析程序，利用Sierpinski地毯进行标定，检验了程序的可行性和准确性。
(2)分割阈值的选取对分形维数的计算具有重要影响，随着阈值的增大，分形维数总体上呈逐渐增大的趋势；给出的分割阈值选取的方法，能够使密封环端面图像中的孔隙结构最大化地表现出来，提高了分形维数计算的精确度。
(3)密封端面形貌具有高度的自相似特征，利用密封端面共聚焦扫描图像结合分形理论，可以实现对密封端面形貌的快速、准确的定量描述，这对于揭示密封端面的泄漏通道具有重要的意义。
