煤是一种固体可燃有机岩。宏观上，烟煤宏观煤岩类型划分为光亮煤、半亮煤、半暗煤和暗淡煤;微观上，煤由镜质组、惰质组、壳质组三大显微组分组及少量矿物组成的[1-4]。煤的化学组成和工业利用中表现的结焦性等工艺性质主要受煤化程度和煤岩显微组分组成影响[5-6]。传统的人工肉眼在显微镜下测量镜质组反射率和鉴定煤岩显微组分，不仅劳动强度大、费时长，而且受人的主观认识差异影响，导致不同操作者测量结果一致性不好，实验室之间的鉴定数据比对难，严重制约了煤岩参数在煤炭分类和炼焦配煤等工业中的应用，实现煤岩自动化测试是多年来煤岩学界迫切盼望解决的技术难题。针对这一问题，国内外开展煤岩自动测试技术研究已有50 余年历史，研究成果主要体现在煤岩显微光度计自动测定技术、煤岩显微图像自动识别技术两个方面，目的是获得煤岩显微组分组含量和镜质体平均反射率[7-22]。但是，在油浸反光条件下，由于粉煤光片中黏结剂与壳质组都呈现低反射特征，存在两者灰度重叠区;当抛光面上黏结剂区域受其下方煤颗粒反射的影响时，抛光面上显示不同亮度、形态模糊，不呈现黏结剂的低反射率和低灰度特征;壳质组与镜质组、镜质组与惰质组之间反射率或灰度存在过渡现象。鉴于这些原因，目前仍然不能很好地实现三大显微组分组的煤岩自动化准确识别测定。为此，煤炭科学研究总院西安研究院立项开展攻关研究，研发成功基于K均值聚类算法的煤岩显微组分图像自动识别和镜质体反射率自动检测系统。需要说明的是，本检测系统适用于烟煤的三大显微组分组识别测定，当镜质体反射率>1.30%时，可识别测定镜质组和惰质组;煤中矿物的自动识别将另行处理。
1 系统结构和主要模块
煤岩显微组分图像自动识别分析系统的开发，涉及煤岩、机械自动化、光电、计算机、图像处理等多学科融合的技术，由显微图像自动采集平台和显微组分组图像自动识别软件子系统组成。
1.1 煤岩显微图像自动采集平台
煤岩显微图像自动采集平台，即煤岩显微图像自动采集子系统，主要包括反射偏光显微镜、工业相机、自动扫描台、自动聚焦装置、相关衔接接口装置和计算机控制软件模块(图1)。计算机控制模块主要包括自动扫描控制子模块、自动聚焦控制子模块及图像采集控制子模块。
图1 煤岩显微图像自动采集平台
Fig.1 Automatic collection platform for coal micro-image
图像自动采集平台运行流程为:
(1)在油浸反光条件下，按设定样本点距及行距，启动自动扫描控制子模块，驱动自动扫描台移动至样品的初始测点;
(2)图像采集控制子模块发出指令，工业相机采集初始测点的第1次显微图像;
(3)将采集的初始测点第1次图像输入自动聚焦控制子模块，判别成像清晰度，发出调焦指令，显微镜Z轴自动聚焦装置驱动物镜调节焦距;
(4)工业相机进行初始测点第2次图像采集，自动聚焦子模块对图像成像质量进行第2次成像清晰度判别和调焦，直至图像成像质量符合要求，输出并存储初始测点图像，完成初始测点的聚焦成像;
(5)驱动自动扫描台移动至第2个测点，进行该测点的自动聚焦成像;
(6)重复上述步骤，直至采集完成整个样片全部测点的图像，结束样品图像的自动采集。
1.2 煤岩显微组分组图像自动识别软件子系统
该子系统包括显微图像读入、显微图像去噪预处理、显微组分组图像自动识别、显微组分组含量分析、镜质体反射率自动检测、检测报告生成输出等模块(图2)。其中，镜质体反射率自动检测模块包括均质镜质体图像自动识别子模块和标准物质灰度-反射率关系模型子模块。
图2 显微组分组自动识别分析子系统结构
Fig.2 Automatic recognition and analysis subsystem for maceral group image
(1)显微图像去噪预处理模块。该模块功能是剔除图像中的黏结剂区域以及抛光面下方受煤颗粒反射影响的黏结剂区域，输出仅包含煤颗粒的有效图像，难点是开发出将煤颗粒显微图像从黏结剂中有效分割出来的显微图像去噪预处理技术。
(2)显微组分组图像自动识别模块。显微组分组图像自动识别技术是煤岩自动化检测的核心技术。运用图像处理技术，实现镜质组(V)、惰质组(I)和壳质组(E)三大显微组分组图像的自动分割、识别。分别输出分割后各显微组分组的图像信息，为下一步定量分析提供数据。
(3)显微组分组含量分析模块。煤是由三大显微组分组及少量矿物组成的，笔者主要针对显微组分组进行图像识别分析技术开发。鉴于矿物的图像灰度特征横跨显微组分组的各灰阶分布范围，其人工鉴定是在干物镜下进行的，故此对其图像的自动识别技术将另行开展研究。显微组分组采用面积统计代替中国国家标准标中的点统计。具体的说，就是先分别统计每一幅有效显微图像中各显微组分组的有效像素和并逐一相加，分别获得样品镜质组、惰质组、壳质组的像素和，同时将3者累加获得样品的总像素和;然后分别用前3者除以后者，分别得到镜质组、惰质组和壳质组的百分含量。统计公式为
(1)
(2)
式中，∑C为各显微组分组的总有效像素数;∑V，∑I，∑E分别为镜质组、惰质组、壳质组的总像素;k为第1幅至第n幅有效测点显微图像。
(4)镜质体反射率自动检测模块。在已输出的显微组分分割图像基础上进行均质镜质体亚组分的图像自动识别，并建立标准物质的灰度与反射率线性关系模型，将均质镜质体图像各像素点灰度代入建立的关系模型，换算出各像素点反射率。
(5)报告处理输出模块。用于完成对显微组分组体积含量、镜质体反射率测试的分析报告文件，包括报告信息的录入、对检测数据、图像、报告的输出、保存、打印等。
2 煤岩显微组分组自动识别关键技术
2.1 煤岩显微图像去噪预处理技术
煤岩显微图像去噪预处理的目的是去除影响显微组分组识别的干扰因素，提取有效的煤颗粒图像区域。主要解决2个问题:一是分割灰度重叠区的黏结剂与壳质组;二是剔除受到抛光面下方煤颗粒反射影响的模糊状、半透明黏结剂图像区域。
针对壳质组与黏结剂分割问题，LESTER Edward[7]运用阈值法及颗粒孔洞充填图像处理技术分割黏结剂与壳质组。由于壳质组与黏结剂存在灰度重叠，仅用阈值法无法将2者有效分割，并且壳质组并非都分布在颗粒内部，对于壳质组与黏结剂直接接触部位用孔洞充填处理无法将壳质组分割开。LESTER Edward[21]尝试对同一测点图像进行多次曝光采集，根据不同曝光时间下壳质组与黏结剂图像变化特征差异进行分割识别。CLOKE Michael等[22]尝试运用壳质组与其他成分的荧光图像特征差异进行分割。这些尝试耗时长、过程复杂、效率低，难以在实际测试工作尤其是工业化检测工作中推广。
针对煤颗粒与黏结剂灰度、形态、纹理等特征，采用多种图像综合处理技术，开发出煤颗粒有效提取的煤岩显微图像去噪预处理技术。采集的原始图像为灰度图像，单幅图像像素为1 392×1 040，图像灰度8 bit(灰阶0～255)。
图像处理算法流程如图3所示。
图3 煤岩显微图像去噪预处理算法流程
Fig.3 Procedure of preprocessing algorithms for denoise of coal micro-image
(1)采用多峰值阈值、OTSU算法对图像进行二值化处理，然后辅之形态学处理方法中的区域形状阈值，判断区域位置边缘，得到初始预处理图像。
(2)根据煤颗粒与黏结剂的形状特征，采用方形度、轮廓不平度、轮廓峭度和离心率等形态指标，以及傅里叶描述子、小波描述子、Hu矩描述子、Zernike矩描述子和小波矩描述子等数学运算形状描述子，描述煤颗粒与黏结剂形状特征。
(3)结合纹理特征量，依据主成分分析方法，提取煤颗粒边界轮廓形状描述的最佳指标体系。
(4)将其最佳指标体系作为特征向量，输入到支持向量机神经网络，进行大量的样本训练学习，完成待测样本的去噪预处理。
(5)将黏结剂及受到抛光面下方煤颗粒反射影响的黏结剂区域与煤颗粒区域分割开，灰度置为0，输出有效煤颗粒图像。
图4为两幅典型较难分割处理的粉煤光片原始显微图像人工标定与运用上述算法进行自动去噪处理效果对比图。图4(a)和(e)为采集原始显微图像，图像4(a)处理难点为图像中壳质组与黏结剂像素相邻且灰度重叠，图像4(e)处理难点为图像中模糊状、半透明黏结剂图像区域的剔除;图4(b)，(f)为人工标定去背景噪声图像，图4(c)，(g)为运用本次开发算法自动去噪处理效果图。从生成的图4(d)和(h)灰度直方图也可看出，自动去噪处理后图像灰度分布与人工标定图像灰度分布非常一致，结果极为相符。
图4 粉煤光片典型背景噪声图像去噪处理效果对比
Fig.4 Image of denoised preprocessing effect with typical background noise of polished grain mounts
2.2 显微组分组图像假边界处理技术
各种显微组分之间的亮度及突起差异，造成各显微组分成像时在边界处普遍存在灰度过渡区，形成假边界。假边界区像素灰度特征随不同显微组分组合而异(图5(a),(b))，如壳质组与惰质组边界区的灰度会出现镜质组的灰度，这显然是不合理的(图5(c),(d)，(e))。因此，假边界区内的像素点灰度不能反映显微组分的真实灰度。
形成显微组分假边界区的组合类型包括:镜质组与壳质组边界、镜质组与惰质组边界、惰质组与壳质组边界、惰质组亚组分边界、镜质组亚组分边界、显微组分与黏结剂边界及显微组分与矿物边界。组分间像素灰度差异越大、突起差异越大，边界过渡带越宽，即假边界区像素点个数越多;反之，边界过渡带则窄，假边界区像素点个数少，同一显微亚组分内部一般无明显边界。
壳质组与惰质组、惰质组与黏结剂假边界区最宽，可达5～8个像素点;壳质组与镜质组、镜质组与惰质组、镜质组与黏结剂假边界区次之，一般4～6个像素点;壳质体与半丝质体、壳质组与黏结剂边界一般3～5个像素点;镜质体与半丝质体边界较窄，一般2～4个像素点;同类显微组分亚组分间分界最窄，多不显示边界。
显微组分组假边界区像素点灰度影响显微组分组分割识别的准确性，应作为背景噪声予以剔除。边缘检测是基于灰度突变来分割图像的最常用方法，显微组分组边界普遍适用于边缘检测中的斜坡模型[23](图5(d))。运用适合的边缘检测技术将显微组分边界区像素分割识别出来，然后将分割后的边界区像素作为背景噪声置0处理。
笔者采用Prewitt算子大小为3×3模板分割识别样本显微组分假边界区像素。总体上，样品中组分假边界区像素占比整体煤颗粒图像的10%～18%，对显微组分组准确识别影响极大。假边界区像素占比与样品中各显微组分的含量及其组合关系相关，一般情况下，低变质阶段含镜质组、惰质组、壳质组3种显微组分组合的样本假边界区像素占比高，含其中两种显微组分组合的样本，或者壳质组、惰质组含量极低的样本假边界区像素占比相对较低;中高变质阶段含镜质组、惰质组两种显微组分组合的样本假边界区像素占比低。
2.3 基于K均质聚类图像分割算法的显微组分组图像自动识别技术
显微组分组的自动识别，难点在于目前尚无法建立不同显微组分组之间依据反射率或灰度大小进行判识的量化标准或阈值。PRATT K C[24]利用样品显微组分所有显微图像像素点灰度数据，绘制出灰度频率分布曲线，应用PeakfitTM程序对该曲线进行拟合和分解，分离出一系列呈正态分布的分段频谱曲线(图6)，以代表过渡物质频谱曲线的峰值作为镜质组与惰质组及镜质组与壳质组边界分割阈值，获得镜质组、惰质组和壳质组三大组分组的百分含量。
图5 显微组分假边界区像素灰度特征
Fig.5 Gray level characteristics of pixels in false boundary region of macerals
图6 显微组分图像灰度直方图峰值拟合
Fig.6 Peak fitting of gray histogram of maceral image
David E PEARSON[25]和O’BRIEN G等[14]将样品所有显微组分图像像素点灰度转换为反射率值，绘制显微组分组和矿物反射率累积频率曲线图(图7)，然后绘出中间近似直线的镜质组段切线2，依据镜质组段两端的拐点信息，分别绘出镜质组段两侧壳质组段和惰质组段的切线1和切线3;再由切线1、切线3与切线2的两个交叉点求得镜质组、惰质组、壳质组+矿物三大部分的百分含量。
图7 显微组分灰度累积频率曲线
Fig.7 Gray-level cumulative frequency curve of macerals
但是，实际工作中发现，采用灰度直方图峰值拟合法时，由于峰值多项式拟合阶数的不确定性，造成拟合后峰值个数及位置的多解性，由此确定的显微组分组含量误差较大。另外，当壳质组或惰质组含量高时，显微组分组灰度频率分布曲线呈多峰或不典型正态，各组分组段的灰度累积频率曲线不呈近似直线而是呈明显弧形，采用反射率累积频率FMR方法求得的显微组分组含量误差较大。为此，笔者开发出了基于K均质聚类图像分割算法的显微组分组自动识别技术。
基于聚类分析的图像分割方法[26]为将图像空间的元素用对应的特征空间点表示，通过将特征空间的点聚类划分，然后再将他们映射回原图像空间以得到分割结果。K均质聚类分割算法[27]为以图像像素点灰度特征作为相似度特征指标，采用误差平方和作为聚类准则函数，把样本空间m×n个象素点划分为K个簇，使簇内具有较高的相似度，而簇间的相似度较低。为了尽量使聚类能够自然收敛，需要对算法迭代策略进行改进。
在实验中观察到那些远离聚类中心的点会大大改变聚类中心，改进的中心策略是减小那些对聚类中心干扰过大点的影响。为此，在算法中加入了高斯分布函数为这些点分配权值，使远离聚类中心的点在进行计算聚类中心时产生小的影响。
高斯改进K-means聚类实现算法如下:
(1)输入初始种子点
和所有要分类的vec(x)。
(2)计算vec(i,j)到最近的聚类中心
的距离。
(3)
vec(x)被分入最近的类别K中标记为veck(x)。
(3)更新分类权重pk(x)。
(4)
(4)更新所有聚类中心。
(5)
其中,veck(x)表示那些属于的向量。
(5)当
输出分类结果veck(x)和聚类中心
否则重复步骤(2)。
“聚类个数K”及“初始种子点C”是影响K均质聚类算法聚类结果的关键参数，分别对应样本中包含的显微组分组类数以及各显微组分组的灰度分布中心。聚类中心个数K的确定，可以结合煤的变质阶段显微组分变化特征解决。当煤的镜质体反射率Rmax≤1.30%时，按镜质组、惰质组、壳质组3类聚类;由于在煤化过程中，壳质组对温度敏感，其光学特征和化学性质变化比镜质组和惰质组快，当Rmax>1.30%时，壳质组的颜色和突起与镜质组趋于一致，仅可根据其光学各项异性结合原始轮廓识别壳质组[28]，本文开发图像采集平台不适于光学各项异性观测，故按镜质组、惰质组2类聚类。“初始种子点C”根据样品显微图像显微组分灰度频率分布图特征自动分析获取。
算法具体实现流程如图8所示，根据不同变质阶段样本图像显微组分组合特征，自动判识聚类个数及初始种子点，经过N次迭代运算后将聚类结果映射回原图像空间，输出显微组分组的判识结果。
K均质聚类算法的实际运行效果如图9所示。原始图像9(a)中包含镜质组、惰质组、壳质组3种显微组分组，壳质组与镜质组、镜质组与惰质组分间之间存在灰度过渡现象(图9(b))。自动识别效果图9(c)～(f)显示，高斯判断函数有效地控制了过渡像素点的聚类偏移问题，显微组分组自动识别结果与人工鉴定结果一致。
图8 显微组分组图像的K-均质聚类分割识别算法流程
Fig.8 Procedure of K-means algorithms for the image identification of maceral group
图9 K均质聚类分割算法运行效果
Fig.9 Running results of K means algorithms for typical images
3 有效测点确定与精度分析
3.1 有效测点的确定
为了确定样品有效测点数，选取5个样品，采用0.3 cm×0.3 cm等间距测点布满全片，按300,500,800,1 000幅4种规格采集图像，并进行图像预处理，显微组分组自动识别结果及与人工检测数据对比分析结果见表1和图10。
可以看出:4种采集方式检测结果的重复性都满足国家标准GBT 8899—2013的要求[29]，检测结果的极差很小，说明显微组分组检测结果接近。显微组分组体积分数P≤10%的极差不大于0.9%，体积分数10%<P≤30%的极差不大于2.7%，体积分数60%<P≤90%的极差不大于3.1%。
4种采集规格下与人工检测结果极差都小于5%。其中，按300幅采集样本极差≤3%占比87%，500,800,1 000幅采集样本极差≤3%占比都为93%。为与人工检测国标方法(>500有效测点)进行比对，本次所有样本有效图像选取设定为500幅。
3.2 测定结果重复性分析
为了对测定结果的重复性、准确性进行验证，在中煤科工集团西安研究院煤岩实验室保留的20世纪70年代制订“中国煤炭分类”(GB 5751—86)标准[5][30]时采集的全国各地区400多个煤样中，选取45个不同变质阶段的代表性煤样，对其粉煤光片重新抛光、显微图像采集，采用本文提出的K均质聚类图像分割法进行了显微组分组自动识别检测，结果见表2。同时，与国内资深煤岩学专家当年人工鉴定的标准结果进行了比对。目前，由于国内外没有成熟的显微组分组自动识别测定技术，无法进行实验室间的再现性分析。
表1 不同规格有效测点下显微组分组含量自动识别测定与人工鉴定结果对比分析
Table 1 Comparison of automatic analysis and manual test of maceral group content under different measuring points
检测方式有效测点1E/%V/%I/%2E/%V/%I/%3E/%V/%I/%4E/%V/%I/%5E/%V/%I/%人工鉴定>5004.5 79.3 16.1 8.9 68.8 22.3 10.7 60.3 29.0 5.0 70.6 24.4 0.1 79.7 20.2 300幅5.6 79.9 14.4 10.7 70.4 18.9 8.9 63.0 28.1 3.4 74.4 22.2 0 80.5 19.5 500幅5.9 78.7 15.4 11.0 71.2 17.8 8.8 61.7 29.5 3.9 72.1 23.9 0 80.6 19.4 自动测定800幅5.4 78.4 16.1 10.2 70.5 19.3 8.0 63.5 28.4 3.8 71.3 24.9 0 81.3 18.7 1 000幅5.3 79.0 15.7 10.4 70.3 19.3 8.1 63.2 28.7 3.8 71.6 24.6 0 80.9 19.1 平均值5.6 79.0 15.4 10.6 70.6 18.8 8.5 62.9 28.7 3.7 72.4 23.9 0 80.8 19.2 极差0.6 1.5 1.7 0.8 0.9 1.5 0.9 1.8 1.4 0.5 3.1 2.7 0 0.8 0.8
图10 不同规格测点下样本显微组分组含量重复性测定
Fig.10 Repeatability test of maceral group content at different measuring points
选取10个代表性样品进行了显微组分组自动识别测定结果的重复性验证(表3)。每个样品进行2次图像采集，分别采用K均值聚类分割算法进行显微组分组的自动识别检测，计算两次检测值的平均值和极差。按照GBT8899—2013中重复性限精密度要求[29]，如果2次检测值极差超出要求，需对该样品进行第3次采集、测试，若3次测值的极差不大于表1重复性限的1.2倍，即视该样品检测结果符合国标精密度要求。
表2 全国煤分类样品显微组分组自动识别测定数据准确度分析
Table 2 Accuracy analysis of automatic recognition test data of maceral group for samples in Chinese classification of coals %
样品号反射率Romax图像法自动识别测定VIE国家标准人工鉴定VIE两方法测值之差dV I E10.57 70.6 21.3 8.2 68.9 24.8 6.3 -1.70 3.50 -1.90 20.67 59.5 24.7 15.9 55.3 26.3 18.4 -4.20 1.60 2.50 30.68 88.8 5.8 5.4 90.9 4.6 4.4 2.10 -1.20 -1.00 40.75 78.0 15.1 6.9 79.3 16.1 4.5 1.30 1.00 -2.40 50.75 74.2 15.3 10.5 73.2 18.3 8.5 -1.00 3.00 -2.00 60.75 86.6 8.4 5.0 83.1 11.7 5.1 -3.50 3.30 0.10 70.88 68.7 20.2 11.1 68.8 22.3 8.9 0.10 2.10 -2.20 80.88 85.0 11.2 3.8 89.0 9.7 1.3 4.00 -1.50 -2.50 90.90 62.6 28.8 8.5 60.3 29.0 10.7 -2.30 0.20 2.20 100.92 46.5 25.7 27.8 43.7 30.0 26.3 -2.80 4.30 -1.50 110.94 55.6 37.1 7.4 60.8 33.4 5.8 5.20 -3.70 -1.60 120.96 67.6 27.7 4.7 72.6 22.1 5.3 5.00 -5.60 0.60 131.00 56.1 39.3 4.6 54.6 40.6 4.8 -1.50 1.30 0.20 141.02 59.1 35.3 5.6 60.9 33.9 5.3 1.80 -1.40 -0.30 151.02 61.7 26.5 11.8 57.0 30.0 13.0 -4.70 3.50 1.20 161.05 50.8 45.4 3.7 56.0 40.0 3.9 5.20 -5.40 0.20 171.05 70.9 26.3 2.8 70.6 24.4 5.0 -0.30 -1.90 2.20 181.11 55.5 42.0 2.5 58.8 39.5 1.8 3.30 -2.50 -0.70 191.12 59.9 32.9 7.3 63.9 24.6 11.5 4.00 -8.30 4.20 201.14 67.7 30.1 2.2 72.0 26.6 1.4 4.30 -3.50 -0.80 211.17 87.3 12.7 —86.8 13.2 —-0.50 0.50 —221.18 74.8 25.2 —70.9 29.1 —-3.90 3.90 —231.18 82.7 17.3 —83.9 16.1 —1.20 -1.20 —241.23 75.9 24.1 —76.8 23.2 —0.90 -0.90 —251.24 88.5 11.5 —84.5 15.0 —-4.00 3.50 —261.27 60.1 39.9 —59.1 40.5 —-1.00 0.60 —271.32 92.8 7.2 —93.2 6.7 —0.40 -0.50 —281.34 75.1 24.9 —76.2 23.8 —1.10 -1.10 —291.35 77.4 22.6 —75.4 24.6 —-2.00 2.00 —301.36 99.0 1.0 —98.4 1.6 —-0.60 0.60 —311.42 86.2 13.8 —82.4 17.4 —-3.80 3.60 —321.45 68.7 31.3 —68.5 31.3 —-0.20 0 —331.47 74.5 25.5 —78.1 21.6 —3.60 -3.90 —341.53 69.2 30.8 —70.3 29.6 —1.10 -1.20 —351.53 81.8 18.2 —79.7 20.2 —-2.10 2.00 —361.57 90.5 9.5 —90.1 9.9 —-0.40 0.40 —371.62 68.6 31.4 —64.7 35.3 —-3.90 3.90 —381.67 90.1 9.9 —93.1 6.9 —3.00 -3.00 —391.69 78.6 21.4 —80.2 19.8 —1.60 -1.60 —401.74 68.0 32.0 —69.2 30.7 —1.17 -1.27 —411.76 82.8 17.2 —82.2 17.8 —-0.60 0.60 —421.77 97.0 3.0 —96.8 3.2 —-0.20 0.20 —431.79 79.6 20.4 —77.3 22.7 —-2.30 2.30 —441.92 97.3 2.7 —96.1 3.9 —-1.20 1.20 —452.25 94.9 5.1 —93.1 6.9 —-1.80 1.80 —平均值差值d-0.32-0.28-0.02标准差值Sd3.0963.0151.872
注:① 在镜质体平均最大反射率小于1.30%时，采用镜质组、惰质组、壳质组3类聚类算法;大于1.30%时，采用镜质组、惰质组2类聚类算法; ② 显微组分百分含量以去矿物基计。
表3 10个样品显微组分组含量自动识别重复性测定结果
Table 3 Repeatability test data of automatic analysis for 10 samples
样品号显微组分组自动测定1/%自动测定2/%平均值/%极差d/%样品号显微组分组自动测定1/%自动测定2/%平均值/%极差d/%E5.4 5.8 5.6 0.4 E7.4 6.5 7.0 0.9 1V88.8 90.0 89.4 1.2 6V55.6 55.5 55.6 0.1 I5.8 4.2 5.0 1.6I37.1 38.1 37.6 1.0 E6.9 5.9 6.4 1.0 E0 0 0 0 2V78.0 79.4 78.7 1.4 7V81.8 80.0 80.9 1.8 I15.1 14.7 14.9 0.4I18.2 20.0 19.1 1.8E5.0 4.4 4.7 0.6E0 0 0 0 3V86.6 86.6 86.6 08V78.6 78.5 78.6 0.1 I8.4 9.0 8.7 0.6I21.4 21.5 21.5 0.1E11.1 9.0 10.1 2.1E0 0 0 0 4V68.7 71.0 69.9 2.3 9V68.0 70.0 69.0 2.0 I20.2 20.0 20.1 0.2 I32.0 30.0 31.0 2.0E8.5 6.7 7.6 1.8E0 0 0 0 5V62.6 61.5 62.1 1.110V82.8 85.0 83.9 2.2 I28.8 31.8 30.3 3.0 I17.2 15.0 16.1 2.2
10个样品20组有效图像显微组分组自动测定值的极差分布显示:显微组分组含量(P)数据范围为P≤10%,10%<P≤30%,30%<P≤60%,60%<P≤90%时，最大极差分别为1.8%,2.2%,3.0%,2.3%，完全符合国家标准相应范围重复性限分别为2.0%,3.0%,4.0%,4.5%的要求(图11)。按显微组分组统计，镜质组含量极差最大为2.5%，惰质组极差最大为3.0%，壳质组极差最大为2.1%。可见，10个样品的重复测定结果复合要求，重复性好。
图11 10个样品显微组分组自动测定值的极差分布
Fig.11 Range distribution of macerals automatic analysis test data for 10 samples
3.3 准确性分析
(1)自动识别与人工鉴定结果准确性比对
对45个代表性样品的K均值聚类分割法自动识别测定，获得镜质组、惰质组、壳质组共110项数据，与人工识别的比对结果见表2。
三大显微组分组的K均值聚类分割法自动识别测定数据与人工鉴定标准结果呈现出非常好的线性关系(图12)。两种方法获取的45个样品三大显微组分组110对测定结果的极差(d)≤5%占比高达95%，镜质组、惰质组、壳质组极差范围分别为0.1%～5.2%,0～8.3%,0.1%～4.2%，平均为2.3%,2.3%,1.5%(图13)。说明K均值聚类分割法自动测定与人工鉴定标准结果一致性很好，准确性很高。
图12 显微组分组自动测定与人工鉴定数据线性关系
Fig.12 Linear relationship between manual and automatic detection data of maceral group content
图13 45个样品显微组分组自动测定值的极差分布
Fig.13 Range distribution between automatic analysis methods and manual test data for 45 samples
(2)偏倚性分析
按照国家标准GB/T18510—2001[31] 准确度分析方法，利用表2中两种方法45个样品数据，计算出镜质组、惰质组、壳质组的平均差值
差值的标准差Sd、统计量tc。统计量tc计算公式为
(6)
计算出的镜质组、惰质组、壳质组统计量tc分别为0.007,0.071,0.431，从GB/T 18510—2001中查得(n-1)自由度下临界值tt分别为2.014,2.014,2.086，即3种显微组分组的统计量tc<tt，说明K均值聚类图像自动识别法与人工鉴定方法没有偏倚。
(3)不同方法的显微组分组识别结果及分析
选取了15个代表性样品，每个样品按固定间距采集625幅有效图像，对采集的原始图像进行了去除黏结剂、去边界处理后，分别采用本文提出的K均质聚类图像分割法(方法1)、直方图峰值拟合法(方法2)、灰度累积频率图FMR法(方法3)对样本进行显微组分组的显微图像自动识别测定，与权威煤岩专家人工鉴定结果作为比对标准，结果数据见表4。
表4 3种自动识别测定与人工鉴定的显微组分组数据对比
Table 4 Comparison of three automatic recognition methods and manual test data of maceral group
样品序号反射率Romax/%显微组分组组合百分含量/%人工鉴定K-means聚类法FMR法峰值拟合法样品序号反射率Romax/%显微组分组组合百分含量/%人工鉴定K-means聚类法FMR法峰值拟合法V 90.9 88.8 77.6 94.2V 58.855.541.865.610.68 I 4.6 5.8 9.5 3.2 91.11I 39.54252.634.4E 4.4 5.4 12.9 2.6 E 1.82.55.60V 79.3 78.0 59.6 81.1 V 63.959.943.370.520.75 I 16.1 15.1 26.1 14.9 101.12I 24.632.945.828.4E 4.5 6.9 14.3 4.0 E 11.57.310.91.1V 68.8 68.7 52.6 71.8 V 76.875.951.281.730.88 I 22.3 20.2 26.9 19.2 111.23I 23.224.148.818.3E 8.9 11.1 20.5 9.0 E 0000V 89.085.067.388.2V 76.275.159.786.440.88I 9.711.219.49.6121.34I 23.824.940.313.6E 1.33.813.32.2E 0000V 60.362.648.967.2V 79.781.864.187.650.90I 29.028.838.429.5131.53I 20.218.235.912.4E 10.78.512.73.3E 0000V 72.667.650.171.8V 69.265.643.373.460.96I 22.127.738.525.2141.74I 30.734.456.726.6E 5.34.711.43E 0000V 5761.746.267.1V 82.282.865.48871.02I 3026.531.721.5151.76I 17.817.234.612E 1311.822.111.4E 0000V 56.050.835.455.281.05I 40.045.458.244.8E 3.93.76.40
为了比较3种不同显微图像自动识别测定结果的准确性，分别对15个样品计算了由3种不同显微图像自动识别测定方法获得的三大显微组分组百分含量与对应的人工鉴定标准结果的极差(d)，结果如图14、表5所示。
图14 3种方法识别结果与人工鉴定数据极差分布
Fig.14 Range distribution of three recognition methods and manual test data
表5 显微组分组3种识别方法与人工鉴定数据极差统计
Table 5 Range statistics of three automatic recognition methods and manual test data of maceral group
显微组分极差d /%K-means聚类法FMR法峰值拟合法镜质组0.1～5.22.710.8～25.918.30.8～10.24.9惰质组0.2～8.32.61.7～26.014.00.1～10.24.3壳质组0.2～4.21.70.6～12.06.60.1～10.43.1
图14、表5给出了方法1、方法2、方法3的平均极差。由此可见，本文提出的K均质聚类图像自动识别测定出的显微组分组结果与人工鉴定标准测定的结果很接近，比其他两种方法的准确性更好。
4 结 论
(1)研制了适于偏反光显微镜的煤岩显微图像自动采集硬件平台，主要包括自动聚焦模块、自动扫描模块和显微图像自动采集模块3个部分，成像效果好、采集速度快;建立了煤岩显微组分组图像自动识别工作流程。
(2)开发出煤岩显微图像去噪预处理技术，有效解决了壳质组与黏结剂灰度重叠的有效分割问题，消除了抛光面下方煤颗粒反射对黏结剂区域灰度影响。
(3)形成了基于Prewitt算子的煤岩显微组分假边界图像处理技术，剔除了占比10%～18%假边界区像素对显微组分组自动识别结果的影响。
(4)创建了基于K均值聚类的煤岩显微组分组图像自动分割和识别技术，形成了烟煤的煤岩显微组分组图像自动识别系统。
(5)比对表明，本文开发系统的自动识别结果与资深煤岩专家人工鉴定的镜质组、惰质组、壳质组的极差平均值分别为2.3%,2.3%,1.5%，说明本文开发出的煤岩显微组分组图像自动识别系统测定结果重复性好、准确性高，与人工鉴定方法没有偏倚。
[1] 武汉地质学院煤田教研室编.煤田地质学(上册)[M].北京:地质出版社,1979.
[2] 杨起,韩德馨.中国煤田地质学(上册)[M].北京:煤炭工业出版社,1980:22-44.
[3] 张群,潘治贵.烟煤的宏观煤岩分类系统研究[J].煤田地质与勘探,1999,27(2):2-5.
ZHANG Qun,PAN Zhigui.Establishment of macroscopic classification system for bituminous coal in China[J].Coal Geology & Exploration,1999,27(2):2-5.
[4] 龚至从,张秀仪.腐植煤的显微组份划分及命名[J].煤炭学报,1979(2):11-20.
GONG Zhicong,ZHANG Xiuyi.Grouping and nomenclature of microscopic constituents of humic coals[J].Journal of China Coal Society,1979(2):11-20.
[5] 张秀仪,龚至从,门桂珍,等.烟煤的成因——工业分类[J].煤炭学报,1981(2):30-37.
ZHANG Xiuyi,GONG Zhicong,MEN Guizhen,et al.Genetic-technical classification of bituminous coal[J].Journal of China Coal Society,1981(2):30-37.
[6] 叶道敏,肖文钊,罗俊文,等.煤岩配煤和焦炭强度的预测[J].燃料与化工,1998,29(5):233-236.
YE Daomin,XIAO Wenzhao,LUO Junwen,et al.Coal petrology blending prediction of coke strength[J].Feul & Chemical Processes,1998,29(5):233-236.
[7] LESTER E,ALLEN M,CLOKE M,et al.An automated image analysis system for major maceral group analysis in coals[J].Fuel,1994,73(11)1729-1734.
[8] CLOKE M,LESTER E,ALLEN M,et al.Repeatability of maceral analysis using image analysis systems[J].Fuel,1995,74(5):654-658.
[9] O’BRIEN G,FIRTH B,ADAIR B.The application of the coal grain analysis method to coal liberation studies[J].International Journal of Coal Preparation and Utilization,2011,31(2):96-111.
[10] 金奎励,夏俭.煤组分组定量与镜质组反射率测定的自动化测试[J].中国矿业学院学报,1986(1):60-67.
JIN Kuili,XIA Jian.The automated coal petrology in determining maceral group composition and the reflectance of vitrinite[J].Journal of China University of Mining & Technology,1986(1):60-67.
[11] 胡德生,王文韬,刘其真.数字化自动煤岩分析技术的开发[J].钢铁,2005,40(7):17-21.
HU Desheng,WANG Wentao,LIU Qizhen.Development of digital automatically analysis technique for maceral specification[J].Iron and Steel,2005,40(7):17-21.
[12] 王文韬,胡德生,尹文义,等.数字化煤岩分析系统的设计与实现[J].中国图象图形学报,2003,8(7):783-787.
WANG Wentao,HU Desheng,YIN Wenyi,et al.Design and implementaion of digital coal petrography analysis system[J].Journal of Image and Graphics,2003,8(7):783-787.
[13] O’BRIEN G,JENKINS B,OFORI P,et al.Semi-automated petrographic assessment of coal by coal grain analysis[J].Minerals Engineering,2007,20(5):428-434.
[14] O’BRIEN G,JENKINS B,ESTERLE J,et al.Coal characterisation by automated coal petrography[J].Fuel,2003,82(9):1067-1073.
[15] GUPTA S,SHEN F,LEE W J,et al.Improving coke strength prediction using automated coal petrography[J].Fuel,2012,94:368-373.
[16] 姚伯元,许国贤.煤/焦反射率自动测定系统与技术方法[J].煤田地质与勘探,1996,24(4):18-21.
YAO Boyuan,XU Guoxian.The tool and technical method for measuring coal/coke reflectance[J].Coal Geology & Exploration,1996,24(4):18-21.
[17] 姚伯元,许国贤.自动确定煤岩参数的数据处理技术[J].煤田地质与勘探,1996,24(5):18-21.
YAO Boyuan,XU Guoxian.The data processing technique used to automatic determining coal petrologic parameters[J].Coal Geology & Exploration,1996,24(5):18-21.
[18] 陈洪博,白向飞,李振涛,等.图像法测定煤岩组分反射率工作曲线的建立与应用[J].煤炭学报,2014,39(3):562-567.
CHEN Hongbo,BAI Xiangfei,LI Zhentao,et al.Working curve establishing and application of determining maceral reflectance by image analysis method[J].Journal of China Coal Society,2014,39(3):562-567.
[19] 王培珍,殷子睆,王高,等.一种基于PCA与RBF-SVM的煤岩显微组分镜质组分类方法[J].煤炭学报,2017,42(4):977-984.
WANG Peizhen,YIN Zihuan,WANG Gao,et al.A classification method of vitrinite for coal macerals based on the PCA and RBF-SVM[J].Journal of China Coal Society,2017,42(4):977-984.
[20] 王素婷,朱宪坤,吕青.基于RILBP-GLCM算法的煤岩显微组分识别[J].煤炭学报,2017,36(3):142-144.
WANG Suting,ZHU Xiankun,LÜ Qing.Coal rock macerals recognition based on RILBP-GLCM algorithm[J].Journal of China Coal Society,2017,36(3):142-144.
[21] LESTER E,WATTS D,CLOKE M.A novel automated image analysis method for maceral analysis[J].Fuel,2002,81(17):2209-2217.
[22] CLOKE M,LESTER E,ALLEN M,et al.Automated maceral analysis using fluorescence microscopy and image analysis[J].Fuel,1995,74(5):659-669.
[23] GONZALEZ R C,WOODS R E.阮秋琦,译.数字图像处理[M].北京:电子工业出版社,2011:450-467.
[24] PRATT K C.Use of composite and mosaic imaging of polished surfaces to enhance petrographic analysis of coal by image analysis[J].Org.Geochem,1993,20(6):759-768.
[25] PEARSON D E.Probability analysis of blended coking coals[J].International Journal of Coal Geology,1991(1-4),19:109-119.
[26] COLEMAN G B,ANDREWS H C.Image segmentation by clustering[J].Proc IEEE,1979,67(5):773-785.
[27] MACQUEEN J.Some methods for classification and analysis of multivariate observations[J].Proc.5th Berkeley Symp.Math.Statist,1967(1):281-297.
[28] 秦勇.中国高煤级煤的显微岩石学特征及结构演化[M].徐州:中国矿业大学出版社,1994:18-47.
[29] GB/T 8899—2013,煤的显微组分组和矿物测定方法[S].
[30] 陈鹏.中国煤炭性质、分类和利用[M].北京:化学工业出版社,2001:184-190.
[31] GB/T 18510—2001,煤和焦炭试验可替代方法确认准则[S].
