一种基于R型因子分析和概率神经网络的冲击地压危险性等级评价模型

王佳信，周宗红，李克钢，王海泉，付自国，李晓飞

(1.昆明理工大学 国土资源工程学院，昆明 650093；2.中南大学 资源与安全工程学院，长沙 410083)

冲击地压是指存储在井巷或工作面周围岩体中的弹性应变能，由于瞬时释放而产生剧烈破坏的矿山动力现象。在煤矿生产活动中，一旦发生冲击地压，还可能诱发煤尘、瓦斯爆炸、火灾以及水灾等一系列次生灾害[1-2]。随着浅部煤炭资源开采殆尽，国内外矿山相继进入深部开采，冲击地压的矿井数目越来越多，冲击地压灾害问题越来越突出。目前，国内外专家普遍认为，煤矿冲击地压问题是全世界采煤国家面临的共同难题，因此，冲击地压危险性等级评价显得尤为重要，是有效防治冲击地压的重要基础[3-7]。

近年来，冲击地压危险性预测预报方法概括起来主要有：理论计算、现场监测以及区域性早期评价。理论计算中具有代表性是“载荷三带”理论模型[8-9]，这类方法可以对冲击地压危险性程度进行定量分析，但影响各个采场的采动压力众多，目前尚没有固定的计算公式，主要依赖于经验，主观性较强。现场监测主要通过监测数据的变化规律所释放出的前兆信息对冲击地压危险性进行预测预报，如电磁辐射[10]、钻屑法[11]、微震法[12-13]以及地音法[14]等，该类方法所监测数据往往离散型的，数据容易缺失，冲击地压前兆信息严重缺乏，此外，临界指标的确定是难题。区域性早期评价主要指多参量综合评价，如FDA(Fisher Discriminant Analysis)模型[15]、CPM(Catastrophe Progression Method)模型[16]以及BDA(Bayes Discriminant Analysis)模型[17]，该类方法从不同侧面反映影响冲击地压危险性不确定因素，但需要大量冲击地压危险性评价资料。冲击地压作为一种复杂非线性动力学现象，其机理复杂、影响因素多，冲击地压危险性评价还待于进一步研究，正因如此，应寻求多种方法结合辨识冲击地压危险性等级。

本文提出一种R型因子分析和概率神经网络(Probabilistic Neural Network，PNN)的冲击地压危险性等级评价模型；PNN是在径向基神经网络的基础上融合Bayes 决策理论以及密度函数估计，但是PNN径向基函数还保留径向基神经网络中常采用的高斯函数。基于此，文献[18]采用Alpha稳定分布改进PNN样本层中径向基函数；而文献[19]则利用遗传算法获取模式层矢量最佳数目及匹配的平滑因子参数集，以上的改进都取得了一定的成果。本文借签一种多元统计方法：R型因子分析。采用R型因子析法对冲击地压危险性评价指标数据进行压缩和特征信息提取，即用少数变量尽可能多地反映原有的变量信息，保证原信息损失小且变量数目尽可能少，解决了指标信息重叠问题，此外，经R型因子分析处理后的因子得分数据(新样本数据)PNN模型样本层中基函数有更好的表达。以此为基础，利用因子得分数据建立冲击地压危险性等级评价PNN模型，并以重庆砚石台煤矿冲击地压实例验证模型的可行性。此外，因子分析在工程领域内应用较少，其应用主要体现在指标权重的确定上，目前，对于冲击地压危险性等级评价至今还鲜见文献报道；文章采用R型因子分析对指标降维同时，采用R型因子分析对冲击地压危险性等级进行定性评价。

1 R型因子分析和概率神经网络模型

1.1 R型因子分析模型

因子分析是主成分分析的推广和发展，由相关性矩阵依赖关系出发，将具有错综复杂关系的多个变量(或样品)综合为数量较少几个因子的一种多元统计分析方法。常用的因子分析由两种，一种是R型因子分析(对象为变量)，另外一种是Q型因子分析(对象为样品)；本文采用R型因子分析。

设n个原始变量X1，X2，，Xn，n个变量可由m个因子F1，F2，，Fm与一个n×m阶的系数矩阵A的乘积再加上一个特殊因子ε=(ε1，ε2，，εn)表示(n≥m)

(1)

式中：F1，F2，，Fn为m个公因子；εi为变量Xi(i=1，2，，n)所独有的特殊因子，它们都是不可观测的隐变量，称aij(i=1，2，，n；j=1，2，，m，n)为变量Xi在公因子Fi上是载荷，它反映公因子对变量的重要程度，对解释公因子具有重要的作用。

对式(1)作如下假设：

①公因子彼此不相关，且具有单位方差，即E(F)=0m×1，var(F)=Im×m；

③公因子和特殊因子彼此不相关，即cov(F,ε)=0m×n。

为了消除指标产生的不可公度性，在进行因子分析前需对指标进行标准化处理，具体做法为

(2)

限于篇幅，本文不再列出因子分析的具体步骤，其具体步骤见文献[20-21]。

1.2 概率神经网络模型

PNN是一种基于径向基函数和经典的概率密度估计原理而建立的神经网络[22]，其网络结构见图1，PNN的算法步骤为

图1 概率神经网络结构Fig.1 Structure of probabilistic neural network

首先，将待测样本向量X输入输入层，其中神经元数目与样本维数相等。样本层(部分学者称为模式层)计算待测样本向量X与训练样本间的距离，该层每个节点单元的输出计算为

f(X,Wi)=exp[-(X-Wi)T(X-Wi)/2δ2]

(3)

式中：X为待测样本向量；Wi为输入层到样本层的权重；δ为平滑参数。

然后，对求和层进行某类的概率密度函数(Probability Density Function，PDF)求和，由Parzen窗(即 kernel函数)方法可得各类PDF估计

(4)

式中：fa(X)为a类的概率密度函数值；Xai为i个训练样本向量；s训练样本个数；r为待测样本向量X和训练样本向量Xai维数。

最后，竞争层输出各类概率密度函数，概率最大值的那一类为1，其他类别为0。

2 冲击地压危险性等级评价的R型因子分析和概率神经网络模型

2.1 因子分析的概率神经网络模型预测步骤

冲击地压危险性等级评价的R型因子分析和PNN模型技术路线，如图 2 所示，主要评价步骤如下：①采用R型因子分析对指标数据进行降维处理；②将R型因子分析后的因子得分数据作为PNN输入层，建立冲击地压危险性等级评价的PNN模型；③设置径向基函数的扩展系数，采用SPREAD表示；④得出冲击地压危险性等级评价结果。

2.2 建立冲击地压危险性指标体系

根据砚石台煤矿30多年的动力现象资料，并结合周健等、文畅平以及代高飞的研究成果，选取煤厚(X1)、倾角(X2)、埋深(X3)、构造情况(X4)、倾角变化(X5)、煤厚变化(X6)、瓦斯浓度(X7)、顶板管理(X8)、卸压(X9)以及响煤炮声(X10)等 10 个影响因素作为冲击地压的评价指标。其中，X4，X5，X6，X8，X9以及X10等状态参量指标需要进行数量化处理(见表1)，其数量化参考文献[23]，并将砚石台煤矿冲击地压分为微冲击、弱冲击、中等冲击和强冲击4类，分别用G1，G2，G3和G4表示，其学习样本见表2和表3。

图2 R型因子分析的PNN冲击地压评价模型技术路线Fig.2 Technology Roadmap of rock burst evaluation model of PNN based on R-type factor analysis

2.3 R型因子分析结果

采用SPSS统计软件(自动执行指标数据标准化)对表2中35个样本10个指标数据进行因子分析，其中，主因子数目的确定，要求尽量提取多的信息且又能够达到降维目的。在不知道考虑几个公共因子的情况下，一般优先考虑较多的公因子，然后根据结果再减少因子数，但此过程过于繁琐。Mardia给出不同公共主因子数，所应具备最少的原始变量数之间的关系(见表4)；国内学者常采用累积方差贡献率大于85%(部分学者采用70%～80%)来确定主因子数目。文章首先将主因子个数定为5进行因子分析。因子分析的总方差解释见表5。由表5可以看出，经最大方差法旋转以后，前5个因子变量的特征值均大于1，前5个主因子累积方差贡献率为 89.633%，包含原有指标数据信息89.633%，提取前5个主因子较为合适。

表1 冲击地压危险性等级评价及定性指标赋值Tab.1 Classification of rock burst risk evaluation and qualitative index assignment

表2 样本表Tab.2 Training samples

表3 冲击地压发生地点Tab.3 Location of rockburst

主因子载荷表征的是主因子与原始变量之间的相关系数，主因子载荷矩阵旋转之后载荷系数(见表6)更接近1或者0，这样得到的主因子能够更好的解释和命名变量。结合贡献率与正负相关性作用，由表6可以看出，第一主因子F1与指标X4，X8，X9和X10显著正相关，因此，F1综合构造情况、顶板管理、卸压以及响煤炮声等指标信息。第二主因子F2与指标X2和X5显著正相关，F2综合倾角和倾角变化指标信息，F2可称为倾角因子。第三主因子F3与X6和X7变量存在明显正相关性，因此，F3综合了煤厚变化和瓦斯浓度指标信息。第四主因子F4仅与X3显著正相关，F4可称为埋深因子。第五主因子F5仅与X1显著正相关，F5可称为煤厚因子。

表4 主因子数与原始变量数关系Tab.4 The relationship between the number of principal factors and the number of original variables

(5)

(6)

(7)

(8)

(9)

将表2中样本数据标准化，然后，将处理后数据代入式(5)～式(9)进行主因子得分数据的计算(见表8)。根据主因子得分数据。目前，许多专家采用主成分分析得分图反映对象与指标之间关系[24]，因子分析作为主成分分析的推广，同样可采用主因子得分图反映对象与指标之间关系，作出F1与F2间主因子得分散点图(见图3)。

结合由图3和表6可以看出，主因子F1与F2分别包含原来信息量的31.595%和20.676%；编号为16～22，26～28，30，33，和35的13个样本分布在F1和F2第一象限(正向区间)，第一、二主因子得分较大，这13个样本的煤层倾角较大，多以急倾斜煤层为主，倾角变化较大，构造情况较为复杂，顶板管理较差，支护差甚至无支护，卸压效果差，有的甚至无卸压措施，响煤炮声较多，这与表2中指标实际情况大体一致，可以大致分析这13个样本的冲击地压多以中等冲击～强冲击地压为主。第四象限中5，6，7，9，12，14，15，23，31，32和34等11个样本分布在F1的正向区间和F2的负向区间；表明这11个样本构造较为复杂，顶板管理较差，支护差甚至无支护，卸压效果一般，有的甚至无卸压措施，响煤炮声较多，煤层倾角多以倾斜为主，倾角变化较小，有的甚至无变化，这11个样本的冲击地压多以微冲击～中冲击地压为主。10，11和13等3个样本分布在F1的负向区间和F2的正向区间(第二区间)，表明这3个样本的煤层倾角较大，多以急倾斜煤层为主，且倾角变化较大，构造较为简单，顶板管理较好，卸压效果较好，响煤炮声整体偏较少，个别样本的响煤炮声多，这3个样本多以弱冲击～中等冲击地压为主。第三区间包含了1，2，3，4，8，24，25和29等8个样本分布在F1和F2第三区间，且落在第一主因子F1和第二主因子F2的负向区间，这8个样本构造较为简单，顶板支护较好，卸压效果较好，大部分无响煤炮声，响煤炮声整体偏较少，倾角较小，多以缓倾斜为主，倾角变化较小，这8个样本多以微冲击～弱等冲击地压为主。

表5 特征值及其贡献率Tab.5 Eigenvalues and variances

表6 旋转后因子载荷矩阵Tab.6 Factor loading matrix after rotation

表7 因子得分系数矩阵Tab.7 Coefficient matrix of factor scores

图3 主因子F1和F2得分散点图Fig.3 Scatter plot of principal factor scores of F1 and F2

以上分析可以看出，F1和F2得分散点图可以为冲击地压危险性定性分析提供一种很好的思路，F1和F2主要包含原始数据6个指标和52.271%信息。除F1和F2，F1和F3以外，其余主因子间得分散点图难以对冲击地压危险性进行定性分析，只能对反映指标变化趋势，原因在于：其余主因子间主要包含原始数据指标个数和原有信息较少(小于50%)，限于篇幅，本文将不再列出进行具体分析。

采用综合评价函数对35个样本冲击地压危险性等级进行综合评价，其综合评价函数为

(10)

式中：η1，η2，η3，η4和η5分别为主因子F1，F2，F3，F4和F5各自方差贡献率，0.896 33为累积方差贡献率。

将式(5)～式(9)代入式(10)得到

(11)

由表5表8和式(10)可以计算得出综合主因子得分值(见表9)，将综合主因子得分值在等距d=[0.905 8+0.887 8]/4=0.448 4下分为4种类型冲击地压。第一类冲击地压：微冲击地压，其综合主因子得分值取值范围为[-0.887 8，-0.439 4)，该区间包括1～4，6，8，19和25号等8个样本，结合表2可以看出，除6和8号样本被划分到微冲击类外，其余样本的评价结果与实际情况一致。第二类冲击地压：弱冲击地压，其综合主因子得分值取值范围为[0.439 4，-0.009)，该区间包括5，7，10和24号样本，除10号样本被划分到弱冲击类外，5，7和24号样本评价结果与实际情况一致。第三类冲击地压：中等冲击地压，其综合主因子得分值取值范围为[-0.009，0.457 4)，该区间包括9～15，20～23，26～28和30～35号样本，20～23，26，28，30和33号样本的评价结果与实际存在一定偏差外，其余样本的评价与实际情况一致。第四类冲击地压：强冲击地压，其综合主因子得分取值范围[0.457 4，0.905 8]，区间包括16，17，18和19号样本，从表2可以看出，评价结果与实际情况一致。

表9 综合主因子得分及其排名Tab.9 Comprehensive the principal factor scores and their rankings

2.4 PNN模型仿真结果

对表2进行R型因子分析，采用5个主因子得分数据(见表8)作为新的样本数据(包含原始数据信息的89.633%)进行概率神经网络仿真。冲击地压危险性等级评价只是一个理论结果，为了简明描述冲击地压危险性等级，文章将PNN模型竞争层神经元的期望输出值设置为1(G1)，2(G2)，3(G3)和4(G4)。其中，G1，G2，G3和G4代表冲击地压危险性等级。

此外，为检验PNN模型泛化能力，保证冲击地压危险性等级评价结果的可靠性，应避免主观地选取训练样本多、测试样本少的情况。将表7中35个样本构造5种不同的情况进行学习，即将训练样本与测试样本个数比分别设为31∶4，29∶6，27∶8，25∶10和23∶12进行学习，以探讨训练样本个数对PNN模型预测精度影响。其中，SPREAD的值设置为0.5，以保证冲击地压危险性等级得到准确的评价。与此同时，将未经因子分析的样本数据(见表2)同样构造上述5种情况进行学习。

限于篇幅，仅列出训练与测试样本个数比为31∶4的学习情况(图4～图7)。由图4和图5可以看出，指标数据经因子分析降维处理后，PNN模型对1～31号样本进行训练，判别无错误，32～35号测试样本中34和35号样本判别错误，判别结果与实际情况相比，其相对误差均为33.33%，35个样本判别结果与实际结果的平均相对误差为1.90%；经因子分析后的PNN模型其正判率为94.29%。由图6和图7可以计算出，未经因子分析的样本数据(即原始数据)，PNN模型仿真结果中训练样本判别无错误，4个测试样本中，32，34和35号样本判别错误，判别结果与实际情况相比，其相对误差分别为-66.67%，-66.67%和-33.33%，35个样本的判别结果与实际结果的平均相对误差为-4.76%；未经因子分析后的PNN模型其正判率为91.43%。经因子分析后的PNN模型其预测精度同比提高了3.13%。

为进行比较，表2中列出FDA模型、CPM模型以及BDA模型的评价结果。由表2可以看出，经因子分析后的PNN模型，其判别结果与FDA模型、CPM模型以及BDA模型的评价结果基本一致。PNN模型对训练与测试样本个数比分别为31∶4，29∶6，27∶8，25∶10和23∶12的学习统计结果见图8和表10和表11。由图8可以看出，随着训练样本个数的减少，PNN模型的预测精度越来越低；由图8及表10和表11可以看出，经因子分析降维后的PNN 模型，比采用原始指标数据建立的PNN模型更优；经因子分析PNN模型评价结果与实际冲击地压发生情况对照如表 2所示。从表2可以看出，经因子分析PNN模型评价结果与实际情况较为一致，说明因子分析后PNN模型在冲击地压危险性等级评价中是可行的。

图4 经因子分析冲击地压样本1#～31#训练和误差结果Fig.4 Training effects and the error results of rock burst samples 1#-31# by factor analysis

图5 经因子分析冲击地压样本32#～35#测试效果Fig. 5 Test results of rock burst samples 32#-35# by factor analys1is

图6 未经因子分析冲击地压样本1#～31#训练和误差结果Fig.6 Training effects and the error results of rock burst samples 1#-31# without factor analysis

图7 未经因子分析冲击地压样本32#～35#测试效果Fig. 7 Test results of rock burst samples 32#-35# without factor analys1is

图8 训练样本数对PNN模型的影响Fig.8 Effect of training samples on PNN model

3 讨 论

(1)表2中35组样本数据适合采用因子分析？

因子分析作为主成分析的推广，同样可采用KMO(Kaiser Meyer Olkin)统计量和Bartlett P球形检验作为判别因子分析使用条件(见表12)。采用SPSS软件得到KMO和 Bartlett 检验结果(见表13)。由表12和表13可以看出，本文采用因子分析对表2中35个样本10个指标数据进行降维处理是可行的。

(2)SPREAD的参数值影响着PNN模型的预测精度，SPREAD值越小，对函数的逼近越精确，SPREAD值越大，模型预测误差越大，其默认值为1。针对表10和表11的评价结果，本文将PNN模型参数值设置为0.5，使得PNN模型评价效果达到最优。以表8(主因子得分数据)中35个样本数据进行学习，采用1～31 号样本作训练，32～35样本作测试，将SPREAD值设为1.5，以检验PNN模型的可靠性(图9和图10)；由图9和图10可以看出，当SPREAD值为1.5时，模型正判率为91.43%，可见，经因子分析后的PNN模型在冲击地压危险性等级评价中是可行的。

表10 经因子分析后PNN模型的5种学习情况Tab.10 Five learning situations of PNN model by factor analysis

表11 未经因子分析PNN模型的5种学习情况Tab.11 Five learning situations of PNN model without factor analysis

表12 KMO 和 Bartlett 的检验标准Tab.12 Test standard of KMO and Bartlett

表13 KMO 和 Bartlett 的检验结果Tab.13 Result of KMO and Bartlett test

图9 SPREAD值为1.5下PNN模型训练和误差结果Fig.9 Training effects and the error results when the value of SPREAD is 1.5

(3)从R型因子分析-PNN耦合模型仿真结果来看，一些样本的评价结果与实际存在一定的偏差，可能的原因是：①样本容量有限，指标的选取以及量化还不够合理，使得R型因子分析-PNN耦合模型的评价结果与实际存在一定的偏差。②根据钱鸣高等的研究，抛出煤量在10～50 t，震级在1～2级的冲击地压定义为中等冲击；抛出煤量在50 t以上，震级在2级以上的冲击地压定义为强等冲击。从表2表3可以看出，34*号样本位于3 605 N 段切割眼，喷出煤岩50 t，理论上来看，冲击地压类型定义为为强冲击地压。但实际工程中，针对喷出煤岩吨数的测定，可能因方法自身的局限性，研究人员在对测量值取舍时，犹豫再三，往往从安全的角度出发，导致冲击地压等级评价偏高。此外，35*号样本评价过高，但从工程安全角度来说，评价结果是可以接受的。

图10 SPREAD值为1.5下PNN模型测试效果Fig.10 Test results of PNN model when the value of SPREAD is 1.5

(4)不同的变量尺度会对因子分析评价结果产生影响，因为变量尺度的变化会改变变量的方差和均值。为了消除变量产生的不可公度性，因子分析中采用标准化、均值化或极差变化等方法对变量进行无量纲化处理，实际上，所有无量纲化都会导致指标数据信息的损失，这可通过数学得到证明。以均值化为例，设原始数据矩阵X=(Xij)n×p，无量纲化得到矩阵Z=(Zij)n×p，

其中，均值化公式为

(12)

(13)

(5)因子分析后的指标数据(因子得分数据)，PNN模型中高斯函数有更好的表达，即因子得分数据更多服从正态分布(正态分布的概率密度函数是高斯函数的一个实例)。

文章对原始指标数据和因子得分数据进行正态性分析，首先提出统计假设：H0为指标数据服从正态分布；H1为指标数据不服从正态分布。

采用SPSS软件计算得到正态性检验结果，如表14和表15所示，概率P>0.05时，不能拒绝原检验假设H0，认为指标数据服从正态分布。由表14和表15可以看出，原始数据中X2和X3指标数据服从正态分布，假设每组指标数据携带的信息均等，原始数据中有20%的数据服从正态分布；若不考虑因子分析后指标信息的损失，因子得分数据中有60%数据数据服从正态分布，也就是说，因子分析后指标数据高斯函数有更好的数学表达。

值得一提的是，冲击地压不仅受地质条件和开采技术条件的影响，而且还受Ⅱ 类开采技术条件(例如爆破作业、巷帮扩修、推采速率以及工作面停复产等)的影响[26]，使得冲击地压机理及危险性等级评价成为一项艰难课题，冲击地压危险性评价还有待于进一步研究和探讨。

表14 原始数据正态性检验结果(P值)Tab.14 Normality test results of original data(P value)

表15 主因子得分数据正态性检验的P值结果Tab.15 Normality test results of principal factor scores(P value)

4 结 论

(1)冲击地压危险性等级评价模型——R型因子分析和概率网络模型，PNN过程简单，总收敛于Bayes优化解，稳定性高，不需要确定模型的隐含层，样本追加能力强，训练样本可直接赋值给网络，训练时间仅略大于数据读取的时间，不存在局部最优值，可容忍一些判别错误的样本。

(2)采用R型因子分析对冲击地压10个指标进行降维处理，提取5个新的综合指标，尽可能多地保留原始变量的信息，而且提取的综合指标彼此独立，解决了指标信息重叠以及多重共线性问题，此外，降维后的指标数据，PNN模型中高斯函数有更广泛的数学表达，进一步提高PNN模型预测精度。

(3)重庆砚石台煤矿冲击地压危险性等级评价结果表明：指标数据经R型因子分析降维后，5种不同学习情况下PNN模型仍具有较好的评价效果，其正判率分别为94.29%，88.57%，88.57%，85.71%和82.86%；说明R型因子分析和PNN模型结合在冲击地压危险性等级评价中是可行的。