基于主成分回归的自然发火程度判识模型

杜 斌

（1.煤炭科学研究总院，北京100013；2.中煤科工集团沈阳研究院有限公司，辽宁 抚顺113122；3.煤矿安全技术国家重点实验室，辽宁 抚顺113122）

采空区遗煤自然发火是影响我国矿井安全生产的重大危险源之一。当工作面上隅角及回风流中的CO 体积分数出现上升趋势时，需要对自然发火程度进行辨识及预测，以便及时采取有效应对措施；同样，封闭火区启封前也需要对火区内遗煤的自然发火程度进行判识。目前主要是建立自然发火标志性气体预测预报指标体系[1]，以CO 及其派生指标、C2H4、C2H2为主指标，以链烷比和烯烷比以及温度为辅助指标，通过气体分析法对煤自然发火程度进行判识，贾宝山[2]和董绍朴[3]等通过程序升温实验和主成分分析的方法对自然发火标志性气体进行了优选。这些研究成果都表明了自然发火标志性气体可以用于判识煤的氧化程度。林棉金[4]研究了自然发火过程中生成气体的分阶段宏观特性，杨波等[5-6]研究了不同氧化程度煤在氧化蓄热升温过程中的细观孔隙-裂隙结构演变规律，王彩萍等[7-8]研究了煤氧化过程中微观分子结构和官能团的变化规律。这些研究都指出了煤氧化过程中的分阶段特性，但是对划分各阶段的临界点及划分方法尚未提出。尤其是目前尚未提出煤氧化程度的定量化评价方法。为此通过对煤样进行程序升温实验，借助于多元统计方法和SPSS 统计分析软件，对实验数据进行相关性分析、主成分分析及回归，构建自然发火程度预测模型。

1 煤样程序升温实验

1.1 煤样的制备与实验设计

煤样采用变质程度较低的褐煤，采集自麦垛山煤矿6#煤层工作面。按照国标[9-10]的要求，对煤样进行粉碎、研磨、筛分，筛选出粒径小于100 目（0.15 mm）的煤样封存于密封袋中。对煤样进行工业分析，煤样工业分析结果见表1。

表1 煤样工业分析结果Table 1 Industrial analysis results of the coal sample%

实验过程为：①取1 g 制备好的煤样，放入程序升温实验装置的煤样罐中，连接好进气和出气管路；②升温速率设置：25～80 ℃为0.5 ℃/min，80～200 ℃为1.0 ℃/min，200～300 ℃为2.0 ℃/min；③供气流量100 mL/min，用气相色谱仪分析实验生成的气体组分及体积分数，每隔20 min 记录1 次实验结果。

1.2 实验结果

程序升温实验气体产物分析结果见表2。

表2 程序升温实验气体产物分析结果Table 2 Analysis results of gases generated in the programmed heating experiment

从实验结果可以发现，从448 ℃开始，O2体积分数出现断崖式下降，说明此时已经进入燃烧期，煤样明火燃烧时消耗了大量O2；同时煤体温度先上升至531 ℃后逐渐下降，说明煤样明火燃烧后进入阴燃阶段，此时耗氧量逐渐减小，氧体积分数逐渐升高。待O2体积分数上升回复至与空气中O2体积分数相同时，说明煤样完全熄灭。

N2体积分数在实验过程中基本无变化，这主要是因为N2主要来自于外部环境，其在实验过程中既不是反应物也不是生成物，所以N2体积分数主要是受气体体积相对比的影响而上下浮动。

除了O2、N2之外，CxHy和CxOy的变化趋势大致相同，具有较强的相关性。

2 相关性分析

目前对自然发火程度进行预测主要是通过气体分析法和测温法，由于测温法主要是借助温度传感器等，无法进行非接触式测量，因此气体分析法使用较为广泛。

但是在程序升温过程中产生的其他CxHy和CxOy气体与温度的相关程度大小尚未被定量描述，因此可以将各种气体看作随机向量，通过相关性分析，筛选出与温度相关系数大的气体作为描述自然发火进程的标志性气体。

标志性气体主要考虑从与温度具有较强相关性的CxHy和CxOy中选取，但是CxHy和CxOy的体积分数量级变化范围较大，涵盖了10-1～10-6，需要将实验结果进行标准化变换。标准化变换如下[11]：

计算相关系数矩阵R，由此绘制的相关系数矩阵图如图1。

图1 相关系数矩阵图Fig.1 Correlation coefficient matrix

图1 中蓝色代表负相关，红色代表正相关。圆圈面积的大小表示了相关性的强弱。与温度相关系数最大的是CO，温度与O2、N2呈负相关，与其他气体具有较强的正相关性。

3 主成分分析

CO 在自然发火过程中出现时间最早且贯穿整个自然发火进程，便于检测，因此常常作为预测预报自然发火的标志性气体。但是由于自然发火过程中产生的各气体与温度T 都具有较强的相关性，使用单一气体难以全面刻画自然发火进程，因此需要借助主成分分析，在尽量不丢失描述自然发火进程信息的情况下，将具有相关性的各气体转化为线性无关的主成分。

3.1 主成分分析的基本思想

主成分分析就是设法将原来具有相关性的指标重新组合成一组新的互相无关的几个综合指标来代替原来指标，同时根据实际需要，从中选取几个较少的综合指标尽可能多地反映原来指标的信息[12]。通常数学上的处理就是将原来p 个指标做线性组合，作为新的综合指标。

式中：Xi（i=1，2，…，p）为随机向量，分别代表程序升温实验中出现的气体；Fi（i=1，2，…，p）为主成分；ai1、ai2、…、aip为系数。

求取系数向量（ai1，ai2，…，aip）主要是通过相关系数矩阵求特征值和特征向量，单位特征向量即为系数向量。按照特征值取值大小确定主成分的排名。如果前k 个主成分的贡献率达到了85%，就认为前k个主成分基本包含了全部测量指标所具有的信息。

3.2 主成分分析的结果

借助于SPSS 软件[13-14]的主成分分析功能，计算得到相关系数矩阵的特征值为λ1=6.373，λ2=0.562，λ3=0.049，λ4=0.012，λ5=0.004，λ6=7.837 ×10-5，λ7=3.030×10-5。其中λ1/Σλi=91.038%，因此可以认为第1主成分基本包含了描述自然发火进程的全部信息。第1 主成分的系数向量为特征值λ1对应的单位特征向量，记为a1＝［0.145，0.141，0.154，0.156，0.150，0.151，0.151］。

所以，第1 主成分F1为：

比较系数向量的p 个分量，可以发现a1p都在0.150 附近浮动，说明各气体Xi对第1 主成分的贡献相近。

4 主成分回归

4.1 回归模型的建立

为了对自然发火程度进行预测，还需要将温度T 作为因变量，各气体Xi作为自变量，通过主成分回归分析，构建温度T 和各气体Xi之间的回归预测模型。之前通过主成分分析，已经在尽量保留原来指标信息的情况下将具有相关性的p（p=7）维气体Xi降维，得到了m（m=1）维线性无关的主成分F。这样就简化了回归方程的结构，只需要构建温度T 和主成分F 之间的回归预测模型，然后将主成分F 与各气体Xi的表达式进行回代，就可以得到煤的自然发火程度预测模型。

首先，在XOY 坐标系中绘制散点图，观察温度T 与主成分F 之间的关系，对T 与F 可能符合的回归模型进行初步辨识。主成分F1与因变量T 的散点图如图2。

图2 主成分F1 与因变量T 的散点图Fig.2 Scatter diagram of principal component F1 and dependent variable T

拟合优度R2=0.886，说明拟合效果良好。

然后将式（7）代入上述对数型回归预测模型，可以得到：

式中：X1为CO；X2为CO2；X3为CH4；X4为C2H6；X5为C2H4；X6为C3H8；X7为C3H6。

这样就得到了煤体温度T 和各气体Xi之间的自然发火预测模型。

4.2 回归模型的验证

为了验证对数型模型的有效性，借助SPSS 统计分析软件的曲线估计功能，分别选取线性、二次、三次、S 型、指数型等模型，计算模型的拟合优度R2并和对数型模型的拟合优度R2见表3。

表3 预测模型拟合优度R2Table 3 Goodness of fit（R2）in prediction models

通过对比可以发现，曲线估计功能拟合效果最好的为三次曲线，但其拟合优度依然小于对数型曲线。所以可以将对数型曲线作为自然发火程度的预测模型。

5 结 语

通过开展程序升温实验，分析了自然发火生成的各气体与温度之间的相关性，发现单一气体难以全面刻画自然发火的程度。借助SPSS 统计分析软件，分析了各气体之间的相关性大小，通过主成分分析、主成分回归等方法，发现温度T 和主成分F 之间大致呈对数型函数关系，构建了以温度T 为因变量、气体Xi为自变量的自然发火程度预测数学模型，可以为今后自然发火程度的预测起到一定的借鉴意义。