采空区煤自燃极限参数灰色关联分析及预测

雷昌奎，江莉娟，邓存宝，邓 军，马 砺，王伟峰，张勇敢

（1.太原理工大学 安全与应急管理工程学院，山西 太原 030024；2.陕西山利科技发展有限责任公司，陕西 西安 710075；3.匈牙利佩奇大学 商业与经济学院，匈牙利 佩奇 7626；4.西安科技大学 安全科学与工程学院，陕西 西安 710054）

煤自燃是一个复杂的动态氧化过程，是煤氧复合作用的结果，严重威胁我国矿井的安全生产[1-3]。研究表明，只有当煤体氧化产生的热量大于其向环境散失的热量时，热量的积聚才能导致煤体自燃[4]。煤氧化作用产热与环境散热之间达到平衡状态时的极限外界条件称为煤自燃极限参数，包括最小浮煤厚度、下限氧体积分数和上限漏风强度[5]。由于煤矿开采环境的复杂多变性、煤层地质赋存条件多样性、采空区随着工作面推进动态发展，致使煤自燃极限参数呈现出时变性、非线性、模糊性等特点，因此，采用科学准确的煤自燃极限参数测算方法对采空区煤自燃危险区域判定、高温点范围确定和预防煤自燃灾害的发生有重要意义。

很多专家学者针对煤自燃极限参数测算进行了研究。徐精彩等[6]应用热平衡法推导出煤自燃极限参数的计算方法，并应用BP 神经网络进行了煤自燃极限参数的预测[7]；孟倩等[8]提出了煤自燃极限参数的支持向量机预测模型; 周西华等[9]、郭亚军等[10]以程序升温实验为基础，研究了不同供风量条件下煤自燃极限参数变化；马砺等[11-12]研究了氯盐阻化剂、粒度对煤自燃极限参数的影响；易欣等[13]采用实验方法研究了咪唑类离子液体处理对煤自燃极限参数的影响规律；Ma L 等[14]实验分析了淮南矿区7 个煤矿煤样的煤自燃极限参数；Wang C 等[15]、谭波等[16]研究了不同变质程度煤的自燃极限参数变化规律及特征；Zhang Y 等[17]探讨了煤自然发火升温和降温过程极限参数变化规律；章飞[18]分析了上、下分层采空区遗煤自燃极限参数并划分了煤自燃危险区域；张辛亥等[19]采用分析解法计算煤自燃极限参数，并以凯达煤矿为例进行了验证分析。

综上所述，针对煤自燃极限参数的研究主要集中于实验测试分析，由于影响煤自燃极限参数的因素十分复杂，对这些因素的影响程度以及影响规律等目前了解得还很少，同时，煤自燃极限参数影响因素在时空上呈现动态变化，数学上表现为复杂非线性，在实际计算煤自燃极限参数时通常简化了条件，因而计算结果与现场实际还有一定的差距。为此，在对煤自燃极限参数影响因素分析的基础上，引入灰色综合关联分析方法对煤自燃极限参数及其影响因素之间的关系进行定量分析，找出影响煤自燃极限参数的主要因素；提出一种改进的粒子群优化算法（IPSO）对支持向量回归（SVR）参数进行优化选择，建立参数优化的支持向量机回归模型来预测煤自燃极限参数，从而得到更好的逼近效果，提高预测精度。

1 采空区煤自燃及其影响因素

煤自燃是煤氧化产热和周围环境散热共同作用的结果，是一个复杂的动态热量变化过程。煤自燃的发生受各种复杂因素的影响，是多因素综合作用的结果，除了与自身的氧化放热性相关外，还与所处环境条件有关，采空区煤自燃影响因素如图1。

由图1 可以看出，影响采空区煤自燃的因素十分复杂，从而也决定了影响煤自燃极限参数因素的复杂性，因此，对煤自燃极限参数及其影响因素进行灰色关联分析，确定影响煤自燃极限参数的主要因素和条件，这对煤自然发火的防治十分必要。

图1 采空区煤自燃影响因素Fig.1 Factors of coal spontaneous combustion in goaf

2 采空区煤自燃极限参数灰色综合关联分析

2.1 煤自燃极限参数模型

煤体氧化升温蓄热量大于散热量是保障煤自燃发生的前提，采空区遗煤自燃必须有较高的氧气体积分数，确保遗煤能够发生氧化放热反应；遗煤须达到一定厚度，使产生的热量能够积聚；且漏风强度不能太大，避免产生的热量被风流带走。这些外界条件的极限值就是煤自燃极限参数，其数学模型如下[5]：

上限漏风强度：

最小浮煤厚度：

下限氧气体积分数：

由煤自燃影响因素分析及煤自燃极限参数数学模型可知，影响煤自燃极限参数的因素十分复杂，并且在工作面开采过程中，不同时空条件下的影响因素实时动态变化，为了避免在实际应用时的简化计算带来的误差，十分必要采用机器学习方法挖掘各影响因素与极限参数之间的非线性关系，建立煤自燃极限参数预测模型，以更好指导现场煤自燃防控。

2.2 灰色综合关联分析

在客观世界中，有许多因素之间的关系是灰色的，难易辨别各因素之间关系的密切程度，这样就难以找到主要矛盾和主要特性[20]。灰色因素关联分析，目的是定量地表征诸因素之间的关联程度，从而揭示灰色系统的主要特性。

以文献［7］中表1 和表2 中的25 组数据为基础，以上限漏风强度为系统特征行为序列X0=（x0（1），x0（2），…，x0（n）），以采空区距工作面距离（X1）、氧气体积分数（X2）、煤温（X3）、放热强度（X4）、浮煤厚度（X5）作为相关因素序列Xi=（xi（1），xi（2），…，xi（n））。X0与Xi之间的相关关系即上限漏风强度与影响因素关系如图2。

图2 上限漏风强度与影响因素关系Fig.2 Relationship between the superior intensity limit of air leakage and factors

由图2 可知，上限漏风强度X0与其影响因素Xi之间的线性关系非常弱，无法通过线性分析确定X0与Xi之间的紧密程度。灰色综合关联度分析注重事物发展过程的相似性与相近性，更加客观地定量表达事物之间的关联程度。

灰色综合关联度计算公式为：

式中：ρ0i为X0与Xi的灰色综合关联度，ρ0i∈［0，1］；θ 为灰色综合关联度系数，θ∈［0，1］，一般地，可取θ=0.5；ε0i为X0与Xi的灰色绝对关联度，ε0i∈［0，1］；r0i为X0与Xi的 灰 色 相 对 关 联 度，r0i∈［0，1］。

灰色绝对关联度计算公式为：

其中：

灰色相对关联度计算公式为：

其中，

灰色综合关联度既体现了折线X0与Xi的相似程度，又反映出X0与Xi相对于始点的变化率的接近程度，是较为全面地表征序列之间联系是否精密的1 个定量指标[21]。根据以上理论计算公式，取θ=0.5，得出的各因素序列与上限漏风强度的灰色综合关联度见表1。

由表1 可知，浮煤厚度X5与上限漏风强度的灰色综合关联度最大，氧气体积分数X2与煤温X3次之，放热强度X4与采空区距工作面距离X1最小。说明浮煤厚度对采空区煤自燃的上限漏风强度影响最大，这与预防采空区煤自然发火技术手段之一的提高工作面回采率，减少采空区遗煤与丢煤相吻合。浮煤越薄，越不利于热量的蓄积，煤自热环境更容易被破坏，无法达到煤自然发火的极限条件，从而煤自燃就不会发生。

表1 各因素与上限漏风强度的灰色综合关联度Table 1 The grey synthetic degree between the superior intensity limit of air leakage and factors

3 采空区煤自燃极限参数预测与验证

3.1 支持向量回归机

支持向量机[22]（Support Vector Machine，SVM）是由Vapnik 领导的研究小组于1995 年首先提出的一种基于统计学习理论的VC 维理论和结构风险化原理的机器学习方法，可用于模式分类和非线性回归。SVM 方法是把样本数据“升维”，即将样本空间线性不可分数据映射到线性可分的高维空间中，解决样本空间中的高维度非线性问题[23-24]。低维空间到高维空间映射关系如图3。

图3 低维空间到高维空间映射关系Fig.3 Mapping relationship from low-dimensional space to high-dimensional space

支持向量机方法应用于回归问题，得到支持向量回归机（Support Vector Regression，SVR），在支持向量机回归中，输入样本X 通过非线性映射Φ（X）映射到1 个高维的特征空间，然后在这个特征空间中构造优化的线性回归函数。由于无论样本大小、维数高低，径向基核函数均适用，且收敛域宽，只有1个待优化参数γ，易于实现，为此选用径向基核函数。对于径向基核函数的SVR，惩罚因子C 和核参数γ 是影响SVR 性能的主要参数，SVR 预测结果的准确性和可靠性依赖于（C，γ）的最佳选择[25-27]。为此利用改进的PSO 对SVR 的惩罚因子C 和γ 核参数进行优化选择，来寻求最优的平衡参数对（C，γ）。

3.2 粒子群优化算法改进策略

粒子群优化源于对自然界中鸟类捕食行为的研究，并应用于求解优化问题。在每一次搜索迭代过程中，粒子通过种群个体极值与全局极值更新速度与位置：

式中：Vi、Li分别为第i 个粒子在D 维搜索空间的速度与位置；k 为当前迭代次数；ω 为惯性权重；c1、c2为非负常数，称为加速度因子；r1、r2为分布于（0，1）之间的随机数；Pi、Pg分别为个体极值与全局极值。

粒子群优化算法虽然收敛速度快，但在算法的收敛后期，随着种群多样性的丧失，局部搜索能力不强，易陷入局部最优，为此分别从以下3 个方面进行改进：

1）自适应惯性权重ω。对于适应度函数值优于平均适应度函数值的粒子，其对应的惯性权重较小，使得该粒子局部搜索能力增强，反之对于适应度函数值差于平均适应度函数值的粒子，其对应的惯性权重较大，使得该粒子向全局较好的搜索区域靠拢。因此，根据个体适应度值的不同，惯性权重自适应调整，保证了粒子惯性权重的多样性。

式中：ωmax、ωmin分别为惯性权重最大值0.9、最小值0.4；f 为粒子当前适应度函数值；fmin、favg分别为粒子适应度函数值的最小值、平均值。

2）替代机制。每当全局极值更新之后，将所有粒子按适应适度值大小排序，将最差的一半粒子的位置和速度用最好的一半粒子替代，以增强种群的多样性和粒子的局部搜索能力。

3）网格搜索。为了进一步增强PSO 算法的局部搜索能力，提高算法的收敛精度，在PSO 算法中引入网格搜索。由于网格搜索是一种穷举法，如果将网格搜索置于PSO 内循环，势必造成寻优时间过长[28]。为此，在利用PSO 对SVR 的参数（C，γ）进行优化后，再以PSO 寻优的参数（Cp，γp）为中心，在C 维和γ 维上分别以步长cstep 和rstep 布下网格，然后计算二维网格［0.5Cp，1.5Cp］×［0.5γp，1.5γp］区域内各交点的适应度函数值并与全局极值比较，若优于全局极值，则用该点替换全局最优粒子，否则保持全局最优粒子不变。

步长取值方法：

式中：Cmax、γmax分别为参数C 与γ 取值范围中的最大值。

网格搜索示意图如图4。

图4 网格搜索示意图Fig.4 Schematic diagram of grid search

3.3 基于IPSO-SVR 的煤自燃极限参数预测模型

3.3.1 样本数据

以煤自燃极限参数中的上限漏风强度预测为例，选取影响上限漏风强度的主要因素即煤的放热强度、煤体温度、实测氧气体积分数、采空区距工作面的距离和浮煤厚度作为输入，上限漏风强度作为输出，采用文献［7］表1 中的20 组数据作为训练样本，表2 中的5 组数据作为测试样本。为了取消各维数据间数据级差别，避免因为输入输出数据数量级差别较大而造成预测结果误差较大，对输入输出数据进行归一化处理，归一化公式为：

3.3.2 IPSO-SVR 算法流程

根据SVR 原理和PSO 算法理论，IPSO 优化SVR 参数建立煤自燃极限参数预测模型，IPSO-SVR建模流程如图5。

图5 IPSO-SVR 建模流程Fig.5 Flowchart of IPSO-SVR modelling

IPSO-SVR 算法具体计算步骤如下：

步骤1：训练样本和测试样本的载入及其归一化。

步骤2：粒子群优化算法参数初始化。其中，c1＝1.5，c2＝1.7，C∈［0.01，50］，γ∈［0.01，50］，迭代次数为200，种群规模为20，采用5 折交叉验证。

步骤3：产生初始粒子和速度，以交叉验证均方误差作为适应度函数值。依据适应度值，初始化粒子全局极值与个体极值。

步骤4：依据式（8）计算惯性权重，根据式（7）更新粒子的速度和位置。计算适应度值，更新个体极值与全局极值。

步骤5：执行替代机制。

步骤6：判断是否满足终止条件。若满足，则执行步骤7，否则，返回步骤4。

步骤7：执行网格搜索，输出最优参数对（C，γ）。

步骤8：将最优参数赋给SVR，构建参数优化的SVR 预测模型。

3.4 预测结果

3.4.1 性能指标

为了进一步研究模型的预测性能，在保持训练样本和测试样本不变的条件下，同时利用标准PSOSVR 模型（SPSO-SVR）、多元线性回归模型（MLR）对煤自燃极限参数中的上限漏风强度进行训练建模和测试分析，并与文献［7］、文献［8］的预测结果进行对比。上述模型预测性能评估指标如下：

平均绝对百分误差（MAPE）：

绝对百分误差在3%以内的比例（MP）：

决定系数（R2）：

式中：fi为预测值；yi为真实值；n0为绝对百分误差在3%以内的数量；yˉ为真实值的平均值。

3.4.2 结果对比验证

根据算法设计，经过改进的PSO 对SVR 参数优化选择，寻优结果为C=7.188 4，γ=0.617 8。同时，利用标准PSO 优化SVR，参数寻优结果为C=35.763 2，γ=0.510 1。将寻优参数赋于SVR，建立参数优化的SVR 预测模型，对测试样本进行预测，不同模型预测结果对比结果见表2。

表2 不同模型预测结果对比Table 2 Comparison of prediction results for different models

为了更形象地分析数据，给出几种模型的预测结果，预测结果散点图如图6。

图6 预测结果散点图Fig.6 Scatter plot of the predicted results

由图6 可以看出，IPSO-SVR 预测结果基本与零误差线重合，说明其预测结果与真实值相吻合，具备很好的拟合效果，而其他几种模型相对于离散分布于零误差线周围，尤其是MLR 模型，预测结果距离零误差线较远，说明其预测结果产生了较大误差。不同模型的预测结果相对误差对比如图7。

图7 不同模型的预测结果相对误差对比Fig.7 Comparison of relative error for different models

由图7 可以看出，IPSO-SVR 模型相对误差更小，在2.6%以内，具有更高的预测精度，MLR 模型与文献［7］神经网络模型预测结果明显波动较大，预测结果与真实值之间产生了较大偏差。

不同模型性能指标汇总见表3，不同模型性能对比如图8。

由表3 和图8 可以看出，相比于其他模型，MLR预测结果最差，最大相对误差达到19.788%，MAPE达到9.064 7%，说明线性回归模型用于煤自燃极限参数的预测并不理想；从而也凸显出煤自燃极限参数与其影响因素之间存在更明显的非线性关系，而线性MLR 对非线性关系解释能力有限，不宜用于煤自燃极限参数预测。

表3 不同模型性能指标汇总Table 3 Summary of performance indicators for different models

图8 不同模型性能对比Fig.8 Performance comparison of different models

MAPE＜3%往往作为评估预测模型性能的标准，而工程应用中也以此作为衡量指标[29-30]。总体而言，支持向量机模型的MAPE 均小于3%，而用神经网络模型预测结果的MAPE 为6.15%，同时其MP 和R2均小于支持向量机模型。预测结果表明，与神经网络方法相比，支持向量回归方法预测煤自燃极限参数精度更高、效果更好，更适合于解决小样本、多维非线性映射问题。相比于采用变步长网格搜索方法优化SVR 参数[8]，采用PSO 优化SVR 具有更好的预测效果，而且，IPSO-SVR 预测结果的MAPE 减小达50%以上，说明智能优化算法用于SVR 参数选择具有提高其预测精度的潜力。同时，对标准PSO 进行改进，能够显著提高标准PSO-SVR 模型的非线性处理能力和预测精度。由此可以看出，提出的IPSOSVR 方法能够更加精准地预测煤自燃极限参数。

4 结 语

1）煤自燃极限参数与其影响因素之间的线性关系较弱，为此引入灰色综合关联分析方法对煤自燃极限参数及其影响因素之间的关联程度进行分析，浮煤厚度与上限漏风强度的灰色综合关联度最大，氧气体积分数与煤温次之，放热强度与采空区距工作面距离最小。

2）基于PSO 原理，从迭代初期的自适应惯性权重调整、迭代中期的粒子替代机制和迭代结束后网格再搜索3 个方面对标准PSO 进行了改进，加快PSO 收敛速度，提高其局部搜索能力和稳定性。提出了基于改进PSO 优化SVR 参数的煤自燃极限参数预测模型，并对上限漏风强度进行了预测分析。

3）通过IPSO-SVR、SPSO-SVR、变步长网格搜索优化SVR 参数方法、BP 神经网络和MLR 模型对比分析，IPSO-SVR 可获得较高的预测精度，预测结果的平均相对误差仅为1.140 1%，明显优于其他4种模型。IPSO-SVR 模型具有更精准的预测能力，更好的预测效果，更适合于影响因素复杂的煤自燃极限参数预测。