洪子云， 严承琳， 闵 红， 邢彦军， 李 晨， 刘 曙*

1. 东华大学化学化工与生物工程学院生态纺织教育部重点实验室， 上海 201620 2. 上海海关工业品与原材料检测技术中心， 上海 200135

引 言

煤炭是古代植物埋藏在地下经历了复杂的生物化学和物理化学变化逐渐形成的固体可燃性矿物， 占我国一次能源消费总量的70%以上。 中国是世界上最大的煤炭进口国， 进口煤炭种类包括无烟煤、 烟煤、 褐煤、 其他煤， 在煤炭进口过程中， 为规避国家相关进口管理政策， 套用优惠贸易逃避税收征管， 存在不法企业伪报煤种的现象。 传统煤种鉴别方法需要测定干燥无灰基挥发分、 低煤阶煤透光率、 粘结指数、 恒湿无灰基高位发热量等多项指标[1]， 能耗大， 检测周期长， 不利于口岸快速通关。

煤炭的近红外光谱主要包括N—H， C—H和O—H等含氢基团的倍频与合频吸收， 反映的是煤炭中有机物质及水分的信息， 基于不消耗化学试剂、 快速、 低成本等优势， 采用近红外光谱鉴别煤种受到广泛关注。 有研究针对煤炭近红外谱图， 构建一种置信学习机—支持向量机(CM-SVM)分类器， 对来自山西的肥煤、 焦煤、 瘦煤、 贫瘦煤进行分类， 留一交叉验证准确率为95.48%。 有报道利用深度学习的卷积网络(CNN)提取煤炭样品的近红外光谱特征， 结合极限学习机算法， 测试集最高分类精度可达到96.51%。 宋亮等[2]采集了来自我国主要产煤地的烟煤、 褐煤的可见-近红外光谱， 比较不同算法对于煤种的分类效果， 其中煤体遥感识别(MAO)模型分类时间最短， 随机森林(RF)模型交叉验证判别准确率最高。 Mao等[3]提出了一种改进的多层极限学习机算法， 结合可见-红外光谱， 建立了一种煤种快速鉴别模型， 测试集准确率达到92.25%。 以上研究表明， 近红外光谱特征与煤种具有很好的相关关系， 结合模式识别可用于煤种鉴别， 但目前还未有针对全球不同产地来源煤炭的煤种鉴别应用， 煤炭近红外光谱特征与煤种的相关关系仍有待挖掘。

研究中收集了来自澳大利亚、 俄罗斯、 印度尼西亚、 蒙古、 加拿大、 哈萨克斯坦、 菲律宾、 越南、 朝鲜9个国家410个进口煤炭代表性样品， 涉及烟煤、 褐煤、 无烟煤3个煤种。 对采集的煤炭样品进行漫反射近红外光谱分析， 结合代表性样品成分信息、 X射线衍射分析， 解析了不同煤种煤炭样品的近红外光谱特征及差异。 使用主成分分析(principal component analysis, PCA)对近红外光谱进行数据降维， 结合逐步判别—Fisher判别分析建立煤种鉴别模型， 采用交叉验证和测试样品验证进行模型评价。

1 实验部分

1.1 样品与光谱采集

根据GB 474—2008 《煤样的制备方法》， 将煤炭样品制成粒度小于0.2 mm的化学分析样， 装在密封袋中， 置于阴凉处。 使用Antaris Ⅱ傅里叶变换近红外光谱仪收集光谱数据， 扫描范围为10 000～4 000 cm-1， 波长点数为1 557。 每个样品扫描32次， 平均光谱用作原始光谱。 实验样品来自澳大利亚、 俄罗斯、 印度尼西亚、 蒙古、 加拿大、 哈萨克斯坦、 菲律宾、 越南、 朝鲜这9个国家， 包括烟煤、 褐煤、 无烟煤3个煤种， 共计410个样品， 具体情况如表1所示。

表1 实验样品数量汇总

煤炭样品的漫反射吸光度大小遵循Kubelka-Munk方程， 如式(1)

(1)

式(1)中，KC为吸光系数， 其大小主要取决于煤样的化学组成与分子结构；SC为散射系数， 其大小主要取决于煤样的物理状态， 如颗粒的大小等。

1.2 煤质及XRD分析

为研究煤炭近红外光谱特征与煤种的相关关系， 分别选取2个褐煤、 烟煤、 无烟煤代表性样品， 进行煤质分析及X射线衍射分析。 选取的代表性样品包括印度尼西亚褐煤(YNH)、 菲律宾褐煤(FH)、 澳大利亚烟煤(AY)、 蒙古烟煤(MY)、 俄罗斯无烟煤(EWY)、 朝鲜无烟煤(CWY)。 煤质分析包括近似分析和最终分析， 近似分析根据GB/T 212—2008测定煤炭中水分(Mad)、 灰分(Aad)、 挥发分(Vad)含量， 再根据Mad，Aad，Vad， 计算固定碳(FCad)含量， 计算公式如式(2)

FCad=100-(Mad+Aad+Vad)

(2)

式(2)中：FCad为空气干燥基固定碳的质量分数， %；Mad为一般分析实验煤样水分的质量分数， %；Aad为空气干燥基灰分的质量分数， %；Vad为空气干燥基挥发分的质量分数， %。

最终分析根据GB/T 476—2008， 采用电量重量法测定煤炭中C和H含量， 再计算H/C。 为防止煤炭中的无机矿物质对结构分析产生干扰， 根据GB/T 7560—2001对6个煤炭代表性样品进行脱矿处理， X射线衍射分析实验仪器为德国布鲁克公司生产的D8 Focus X射线衍射仪， 测量条件为： Cu靶辐射， 电压40 kV， 电流 40 mA， 扫描范围为5°～75°， 扫描速度为0.5 s·步-1。 利用PeakFit软件对煤炭样品XRD谱图进行平滑处理， 对位于20°和25°附近的γ峰、 002峰进行拟合， 从拟合曲线中获得峰位(θ)、 面积(α)和半峰全宽(β)等参数。 根据这些参数计算煤炭芳香度(fa)的大小，fa为芳环中碳原子与脂肪族侧链中碳原子的比率， 反映煤炭的芳构化趋势， 计算如式(3)， 式(3)中Aγ和A002分别是对应于γ峰和002峰的峰强度。

fa=car/(car+cal)=A002/(A002+Aγ)

(3)

1.3 建模方法

近红外光谱的变量数通常远远超过样本数， 为解决样本容量小的问题， 建立判别模型前采用PCA进行降维。 PCA由皮尔逊于1901年提出， 霍特林于1933年进一步发展。 PCA主要用于将包含大量非显著成分的多维数据集简化为保留原始数据集中存在的最大变化的少数主成分， 通常用于光谱变量的无监督降维。

逐步判别-Fisher判别分析是一种常用的有监督的分类方法， 先建立已知样本的分类模型， 再对未知样本进行预测。 通过协方差分析的F检验的显著性水平选择最佳的预测变量， 当加入或剔除一个变量后， 对判别函数中的变量进行方差分析， 若计算的F值大于指定的进入值(默认值是3.84)， 则保留该变量； 若计算的F值小于指定的删除值时(默认值为2.71)， 则剔除该变量。 判别的基本原理是将高维度空间的样本投影到低维空间上， 使得投影后的样本数据在新的子空间上有最小的类内距离以及最大的类间距离， 保证最大的可分性。

建模过程中随机选取不同煤种光谱的70%作为训练样本， 剩余30%光谱数据作为测试样本， 如表1所示。 使用判别准确率来评估模型性能， 判别准确率即不同煤种中正确分类的煤炭样本相对于该类样本总数的百分比。

2 结果和讨论

2.1 近红外光谱特征分析

使用Antaris Ⅱ傅里叶变换近红外光谱仪， 采用漫反射积分球模式收集了来自澳大利亚、 俄罗斯、 印度尼西亚、 蒙古、 加拿大、 哈萨克斯坦、 菲律宾、 越南、 朝鲜9个国家， 烟煤、 褐煤、 无烟煤3个煤种， 410个进口煤炭样品的光谱数据， 样品粒度均小于0.2 mm。 图1显示了所有样本在10 000～4 000 cm-1范围内的平均光谱。 总的来说， 不同煤种的近红外谱图差异十分明显， 主要体现在吸光度大小、 光谱斜率、 主要吸收峰三方面。

图1 不同煤种近红外谱图

褐煤的吸光度最低， 大部分褐煤样品的吸光度在0.35～1.1 Abs范围内， 烟煤的吸光度在0.43～1.5 Abs范围内波动， 整体高于褐煤， 无烟煤的吸光度整体最高， 在1.1～1.6 Abs范围内波动； 褐煤的吸光度随波数增加明显降低， 光谱斜率最大， 烟煤次之， 无烟煤的吸光度随波数增加变化缓慢， 光谱斜率几乎为零； 褐煤在5 199 cm-1(O—H伸缩和O—H变形组合振动)处有较强的吸收峰， 在7 021 cm-1(H2O的合频)处有小而明显的吸收峰， 在4 655～4 000 cm-1波段有明显的吸收带(各种有机组合和泛音带)产生[4-5]。 少数烟煤在5 199 cm-1处有较强的水峰， 无烟煤整体谱线趋于平坦， 无明显吸收峰出现。

为了找出造成不同煤种的近红外吸光度大小、 光谱斜率和吸收峰显著差异的原因， 选取褐煤(YNH、 FH)、 烟煤(AY、 MY)、 无烟煤(EWY、 CWY)， 结合煤炭中水分、 挥发分、 固定碳含量、 H/C信息及XRD提供的煤炭芳香度信息进行分析。

2.1.1 吸光度大小

表2是煤质标准分析结果， 展现了不同煤种中水分、 挥发分、 固定碳含量的差异。 图2 (a)是六个代表性样品的近红外谱图， 从图中可以看出不同煤种吸光度大小有显著差异， 无烟煤样品的吸光度最高， 烟煤次之， 褐煤最低。 漫反射吸光度大小主要和KC和SC有关[6]，SC的大小主要取决于样品的物理状态，KC的大小主要取决于样品的化学组成。 实验中所用样品均为粒度小于0.2 mm的粉末， 样本之间物理性质相差不大， 实验中煤炭漫反射吸光度大小主要与KC有关。 而KC大小主要取决于样品种类， 煤化程度越高的煤种， 其芳香环的缩合度越高， 芳香环之间的脂肪侧链不断脱落， 煤的有序结构不断增强[7]， 最终形成完整并具有稳定的共轭环结构和共轭结构的高碳化合物， 导致对近红外的吸收随之增强，KC增大。 图2(b) 是10 000 cm-1处吸光度大小和固定碳含量的关系， 从图2(b)可以看出， 无烟煤固定碳含量最高， 吸光度最大， 褐煤固定碳含量最低， 吸光度最小。 说明煤样近红外漫反射吸光度大小与固定碳含量有关， 固定碳含量越高的煤种， 对近红外光的吸收越强， 整体吸光度越高[8]。

表2 褐煤(YNH、 FH)、 烟煤(AY、 MY)、 无烟煤(EWY、 CWY)样品近似分析和最终分析结果

图2 (a) 褐煤(YNH、 FH)， 烟煤(AY、 MY)， 无烟煤(EWY、 CWY)样品的近红外谱图； (b) 10 000 cm-1处吸光度大小固定碳含量的关系

2.1.2 光谱斜率

六个代表性样品的XRD谱图如图3(a)所示， 图3(a)中每个煤样都具有相同的石墨特征， 并显示出高的背景强度。 表明煤样含有一定比例的无定形碳形式的高度无序的物质[9]。 以EWY样品为例， 图3 (b)显示了EWY在10°～35°范围内位于20°和25°附近的002峰和γ峰的拟合谱图， 通过分峰拟合可以得到峰中心(θ)和半峰全宽(FWHM)等参数， 如表3所示， 根据式(3)计算出芳香度。 采用表2中的H/C反映不同种类样品煤化程度的差异， 无烟煤煤化程度最高， H/C最低。

图3 (a) 褐煤(YNH、 FH)， 烟煤(AY、 MY)， 无烟煤(EWY、 CWY)的XRD谱图； (b) EWY样品在10°～35°范围内的γ峰和002峰的高斯拟合曲线

表3 XRD结构参数

不同煤种fa和H/C的关系如图4(a)所示， 随着H/C的增加，fa在不断降低， 说明煤化程度越高的煤种， 芳构化趋势越明显。 芳构化程度的增加会导致芳香簇和石墨微晶的尺寸不断增大， 煤中的微晶结构逐渐完善， 更接近石墨结构[10]， 芳香分子中的电子跃迁移动到更长的波长， 吸收边缘的产生使得长波方向的吸收系数不断增大[5]， 导致芳构化最强的无烟煤其近红外谱图几乎趋近于水平， 近红外谱图的斜率大小和fa的相关性如图4(b)所示， 可以看出， 随着芳构化程度的增加， 曲线光谱斜率不断减小。

图4 (a) fa和H/C的关系； (b) 光谱斜率和fa的相关性

2.1.3 特征吸收峰

煤炭近红外光谱的主要特征， 是由以脂肪侧链为主的各有机基团对近红外光产生的倍频与合频吸收所导致。 无烟煤由于碳化程度高， 煤样中有机质含量少， 对近红外光吸收少， 无明显吸收特征峰出现， 谱图趋于平坦。 褐煤样品大多为印度尼西亚褐煤， 由于印度尼西亚煤炭大部分来自新生界， 煤层埋藏较浅， 煤炭变质程度普遍较低， 孔裂隙较大， 毛细孔对水的吸附力强， 使得内在水分经空气干燥难以去除[11]， 导致褐煤在7 021 cm-1(H2O的合频)和5 199 cm-1(O—H伸缩和O—H变形组合振动)处有明显的吸收峰， 且5 199 cm-1比7 021 cm-1处的吸收特征更明显。 在5 199 cm-1处产生特征吸收的烟煤样品均为印度尼西亚烟煤， 显示了印啡尼西亚煤炭这个特点。 表2六个样品的水分数据显示， 烟煤和无烟煤的水分含量相差不大， 但褐煤的水分含量明显高于烟煤和无烟煤。 4 655～4 000 cm-1处的吸收峰是有机质光吸收所产生的， 其吸收强弱取决于挥发分大小， 从表2可以看出， 挥发分含量明显褐煤>烟煤>无烟煤， 反映在NIR谱图上吸收峰强弱也是褐煤>烟煤>无烟煤。

2.2 煤种鉴别

上文结合煤炭中水分、 挥发分、 固定碳含量信息及XRD提供的煤炭芳香度信息对煤炭近红外光谱特征与煤种之间的相关性进行了解释， 不同煤种漫反射近红外吸光大小与煤炭中固定碳的含量有关， 近红外光谱斜率与样品芳构化趋势有关， 近红外光谱主要特征吸收峰的强弱与煤炭中水分和挥发分含量有关。 基于这种相关关系， 结合Fisher—逐步判别分析对全球不同产地来源煤炭的煤种鉴别应用进行研究。

首先使用PCA对光谱进行降维， 变量的数量从1 557个减少到394个， 前十个主成分的累计贡献率如图5(a)所示， 其中， PC1 (98.11%)和PC2 (1.71%)这两个主成分解释了99.82%的数据差异， 表明这两个主成分能够最大限度地分离几乎所有样本。 获得的PC1和PC2的载荷图如图5(b)和(c)所示， 载荷图本质上描述的是主成分与全波数变量的相关程度， 载荷系数的绝对值越大， 对主成分的贡献越大。 图5(b)中PC1的载荷图显示， 对PC1影响最大的是位于6 205 cm-1处的特征吸收峰， 其次是位于7 250和5 458 cm-1处的吸收峰。 6 205 cm-1处的特征吸收属于C—H伸缩振动的二倍频吸收， 与煤炭中挥发分对近红外的吸收有关， 7 250 cm-1处的特征吸收属于C—H伸缩振动的二倍频和C—H变形振动的组合吸收峰， 也与挥发分的吸收有关， 5 458 cm-1处的特征吸收属于O—H伸缩振动， 与水分的吸收有关， 说明PC1主要取决于煤炭中挥发分的含量， 其次是水分的含量。 图5 (c)显示PC2与4 345～4 000 cm-1波段的吸光度相关性最强， 其次是位于5 174 cm-1处吸收峰， 4 655～4 000 cm-1波段是各种有机组合和泛音带的吸收带， 与挥发分对近红外光的吸收有关， 5 174 cm-1处的吸收峰表征O—H伸缩和O—H变形组合振动， 与水分对近红外的吸收有关， 说明PC2主要取决于煤炭中挥发分和水分的含量。

图5 (a) 累计贡献率； (b) PC1 (98.11%)的载荷图；(c) PC2 (1.708%)的载荷图

将前十个主成分作为输入变量进行逐步判别， 通过剔除不合适的变量来提高判别准确率。 PC1， PC2， PC3， PC4， PC6， PC7， PC8， PC9和PC10作为有效变量， 进入最终模型。 使用PC1， PC2， PC3， PC4， PC6， PC7， PC8， PC9和PC10建立Fisher判别模型， 共得到2组判别函数和相应的组质心处的函数。 判别函数方程是用非标准系数发展起来的， 要建立一个好的模型， 必须具备高特征值[12]。 判别函数1(57.7%)的特征值为10.63， 判别函数2(42.3%)的特征值为7.802， 均满足条件。 因此， 建立的方程将阐明煤炭样品的良好判别模型。 具体判别函数如式(4)和式(5)所示。

F1=0.301X1+0.212X2-0.211X3-8.169X4+8.917X5-19.393X6-25.721X7-14.089X8+24.106X9-0.079(57.7%)

(4)

F2=-0.031X1+2.001X2+2.083X3-1.830X4-1.799X5-3.920X6+3.020X7+9.314X8-28.885X9-0.003(42.3%)

(5)

式中X1—X9分别代表PC1， PC2， PC3， PC4， PC6， PC7， PC8， PC9和PC10的得分。

将函数1得分作为横坐标， 函数2得分作为纵坐标， 用函数1和函数2的判别得分作散点图， 如图6(a)， 可以看出不同煤种组质心分离明显， 褐煤样品和烟煤样品有部分重叠， 判别函数1对无烟煤和烟煤、 褐煤之间的分离起主要作用， 判别函数2对烟煤和褐煤、 无烟煤之间的分离起主要作用。 判别变量和判变函数之间的相关性可以用绝对相关系数表示[13]， 绝对相关系数越高， 说明变量对判别函数的贡献越大。 图6(b)表明， PC1和第一个判别函数相关性最强， PC2和第二个判别函数相关性最强， PC1和PC2主要与煤炭中挥发分含量有关， 其次是水分， 说明Fisher判别模型对不同煤种进行判别的主要依据是不同煤种中挥发分和水分含量的差异。

图6 (a) 判别函数(F1， F2)得分散点图; (b) 判别函数(F1， F2)中每个变量的相关系数

具体判别结果如表4所示， 建模样品集显示了98%的判别结果， 采用留一交叉验证来验证这些输出， 交叉验证的结果显示97.8%的煤炭样品判别准确， 通过将123个未知样本作为该模型的测试集， 进一步验证了所开发的模型方程， 测试样品集判别准确率为99.1%。 说明近红外光谱结合逐步判别—Fisher判别分析， 能实现对全球范围内煤炭的煤种识别。

表4 判别准确率(%)

3 结 论

利用傅里叶变换近红外光谱收集了来自澳大利亚、 俄罗斯、 印度尼西亚、 蒙古、 加拿大、 哈萨克斯坦、 菲律宾、 越南、 朝鲜9个国家， 涉及烟煤、 褐煤、 无烟煤3个煤种， 共计410个煤种样品的的近红外漫反射光谱数据。 近红外谱图表明， 不同煤种的吸光度、 光谱斜率、 特征吸收峰存在差异。 无烟煤的吸光度整体最高， 烟煤次之， 褐煤最低； 无烟煤光谱斜率几乎趋于零， 而褐煤的光谱斜率最大； 三个煤种中， 仅褐煤的特征吸收峰最明显， 无烟煤几乎无特征吸收峰产生。 为对煤种和近红外光谱特征之间的相关性进行解释， 结合煤炭中水分、 挥发分、 固定碳含量信息及XRD提供的煤炭芳香度信息进行分析， 发现煤炭近红外漫反射吸光度大小与煤炭中固定碳含量呈正相关， 光谱斜率与煤炭芳构化呈负相关， 主要特征峰的强弱与煤炭中水分和挥发分含量呈正相关。 基于近红外光谱特征与煤种之间的相关关系， 结合Fisher—逐步判别分析对全球不同产地来源煤炭的煤种进行鉴别。 首先使用PCA对光谱降维， 然后采用逐步判别筛选出PC1， PC2， PC3， PC4， PC6， PC7， PC8， PC9和PC10作为有效变量建立煤种鉴别的Fisher判别分析模型， 共得到两组判别函数， 判别函数1和判别函数2分别与PC1和PC2相关性最强， PC1和PC2的载荷图又表明PC1和PC2主要与煤炭中挥发分和水分的含量有关， 说明Fisher判别分析对不同煤种进行鉴别的内在依据是煤炭中挥发分和水分含量的差异。 最终测试样品集判别准确率可达到99.5%， 说明近红外结合Fisher—逐步判别分析是一种有效的针对全球范围内煤炭煤种鉴别的方法。