基于DEMATEL-ISM的采空区稳定性模糊评价

Fuzzy Evaluation of Goaf Stability Based on DEMATEL-ISM

摘要

针对采空区稳定性分析影响因素多、定量化评价困难等问题，提出了一种基于决策实验室分析（DEMATEL）—解释结构模型（ISM）的采空区稳定性模糊综合评价方法。利用DEMATEL-ISM模型分析了3个维度共12个主要危险因素间的复杂关系，识别关键因素和因素间的影响路径，进行了因素的层次划分和权重确定。同时，将权重集引入模糊综合评价体系，建立了一种新型采空区稳定性评价模型，实现了采空区稳定性分级。依托某铅锌矿，以930 m中段36-2采空区作为研究对象进行采空区稳定性评价，采用FLAC^3D数值模拟软件分析了该采空区的稳定性状况，并对评价分级结果进行了验证。结果表明：本文提出的采空区稳定性评价模型具有较高的适用性和准确性，能够为采空区的安全管理提供理论支撑。

Abstract

The extraction of mineral resources through techniques such as the room-and-pillar method in underground metal mines results in the formation of numerous goafs，which can readily lead to pillar deformation and roof collapse，thereby creating significant safety hazards.Furthermore，the collapse of goafs can induce mining-related seismic events，water inrush，and mud inrush，posing substantial risks to both personnel and equipment.Conducting precise risk assessments of goaf stability is fundamental for comprehending and managing these risks，which is essential for effective disaster prevention and mitigation.Consequently，performing precise stability risk assessments for goaves is crucial for effective risk management.Addressing the complexities posed by numerous influencing factors and the challenges associated with quantitative assessment in goaf stability analysis，this study introduces a fuzzy comprehensive evaluation approach that integrates the Decision-Making Trial and Evaluation Laboratory （DEMATEL） and Interpretive Structural Modeling （ISM） methodologies.Initially，the factors influencing goaf stability were identified through a comprehensive approach involving field investigations，literature review，and consultations with experts.An initial screening process identified fifteen risk factors associated with rock mass structure，goaf morphology，and environmental impact.Subsequently，twelve critical risk factors were determined，encompassing geological structure，rock mass properties，in-situ stress，groundwater conditions，goaf volume，shape，span，temporal effects，effects of adjacent goafs，mining disturbances，support conditions，and natural disasters.This led to the development of a systematic framework for evaluating goaf stability.Utilizing DEMATEL and ISM analysis，this study investigated the intricate interrelationships among influencing factors to develop a hierarchical framework for goaf stability determinants.An oriented hierarchical structure diagram was constructed to facilitate the classification and weighting of these factors，enabling the calculation of each factor’s significance within the evaluation system.Subsequently，these weights were integrated into a fuzzy comprehensive evaluation system，culminating in the formulation of an innovative model for assessing and grading goaf stability.This model provides management guidelines from an intrinsic safety standpoint，with the objective of fundamentally mitigating the risk of accidents. Through a case study of the 36-2 goaf at the 930-meter level of a lead-zinc mine，the primary influencing factors and direct causes of goaf instability were identified，and the complexity of goaf incidents was thoroughly analyzed.Consequently，the goaf received a stability evaluation rating of “good”. In conclusion，the stability of the goaf was analyzed using the FLAC^3D numerical simulation software.A mine-scale FLAC^3D model was developed to assess the differential stress distribution and failure mechanisms of the rock mass above the goaf，thereby validating the accuracy of the goaf stability classification model.The findings indicate that the goaf stability evaluation model introduced in this study demonstrates high applicability and accuracy.This method mitigates subjectivity in the weight determination process inherent in traditional fuzzy comprehensive evaluation techniques，while simultaneously accounting for a comprehensive range of factors influencing goaf stability.Consequently，this research offers a novel approach to evaluating goaf stability，which may serve as a valuable reference for the stability assessment of goafs in metal mining contexts.

Graphical abstract

关键词

采空区 / 稳定性评价 / DEMATEL-ISM模型 / 模糊综合评价 / 数值模拟 / 稳定性因素

Key words

goaf / stability evaluation / DEMATEL-ISM / fuzzy comprehensive evaluation / numerical simulation / stability factors

随着矿产资源开采规模的不断扩大，采用留矿法和房柱法等空场法开采的地下金属矿山，在开采过程中形成了大量的采空区。采空区的存在极易产生矿柱变形、顶板冒落等安全隐患，不仅会引发空区塌陷事故，而且会诱发矿震、突水和突泥等工程灾害，对作业人员和设备均构成极大的威胁（Wang et al.，2019；解联库等，2020；周科平等，2022；Zhang et al.，2023）。对采空区稳定性风险等级的准确评价是采空区风险认识和控制的基础，对于灾害防控至关重要（Zhang et al.，2021；Wang et al.，2023；蔡鑫等，2024）。

国内外学者提出多种方法对采空区稳定性进行评价，如经验分析、理论分析、数值模拟分析和相似模拟试验等方法（Guo et al.，2019；Jiang et al.，2021；Zhao et al.，2021）。其中，经验分析法凭借实践经验对采空区稳定性进行定性评判，虽然具有一定的实用性，但缺乏精确性和系统性的量化依据。在理论分析法中，蔡泳等（2022）基于Burgers蠕变模型和弹性薄板模型，构建采空区三维蠕变力学模型，推导得到顶板下沉位移量随时间变化的表达式并用于采空区顶板稳定性研究，该方法更侧重于力学计算，对于多种因素作用的采空区顶板稳定与破坏机制的研究认识存在不足。Hu et al.（2012）选用岩石单轴抗压强度、岩石弹性模量等9项指标，基于多元判别分析方法提出了复杂采空区危险程度辨识的贝叶斯判别方法，该方法整合多因素信息，从统计学角度对采空区危险程度进行量化区分，提升了评价的科学性与准确性，但评价过程需要获取大量的现场数据， 且依赖的数学理论模型较为单一，处理高度非线性问题时存在局限性。何标庆（2018）采用数值模拟计算分析采空区群稳定性状况，确定了安全隐患和危险区域；李俊平等（2005）借助相似模拟试验，研究下凹地形采空区顶板自然垮落的最小跨度、移动角和垮落角，揭示了采空区顶板自然垮落特征。模拟手段使抽象的采空区垮落现象得以在实验室环境中呈现并量化研究，但这些方法在模型构建和影响因素选取上仍存在局限性，如难以精准模拟复杂地质条件，影响因素选取不够全面等。

模糊综合评价法是基于模糊数学的评价方法（柯丽华等，2023），用于处理具有模糊性或不确定性因素的评价问题。通过构建模糊集和隶属度函数，将不确定性信息转化为可量化的模糊指标，进而开展综合评价。决策与实验室分析（DEMATEL）是利用图论和矩阵论的原理对系统的影响因素进行分析，确定要素间的因果关系和每个要素在系统中的地位（Dalvi-Esfahani et al.，2019）。郭慧敏等（2023）运用DEMATEL方法分析瓦斯爆炸各致因之间的相互影响关系，并结合ISM对事故致因进行多级递阶层次结构划分，分析导致事故发生的原因要素和结果要素并提出相应的对策措施。解释结构模型（ISM）是一种专为分析系统的相关问题而设计的层级结构分析方法，运用有向图、矩阵、计算机技术和经验知识来分析复杂系统的结构（Warfield et al.，1974）。通过集成DEMATEL和ISM，能够降低ISM中可达矩阵计算的复杂性，减少计算量，使可达矩阵的计算过程变得简单易懂。洪伟斌等（2022）采用DEMATEL和ISM对煤矿透水事故影响因素进行分析，确定了诱发煤矿透水事故的关键因素并提出了对应的管理措施。将DEMATEL-ISM应用于复杂系统中各影响因素之间的相互作用关系分析，能够找出诱发事故的直接影响因素和根本影响因素，并区分影响因素的重要程度。

然而，现有采空区稳定性评价研究仍存在诸多局限性，例如：理论模型和数值模型难以精准模拟复杂地质与多因素之间的耦合影响；影响因素选取不全面，且未能从系统角度剖析因素之间的内在联系，进而影响评价结果的准确性和可靠性；传统模糊综合评价中权重集内影响因素权重划分存在较强的主观性。针对以上问题，采用DEMATEL计算影响采空区稳定性的各因素间的因果关系并确定关键因素；结合ISM 分析层次影响并构建结构图，用于确定影响因素权重，将所得到的影响因素权重应用于模糊综合评价方法，建立基于DEMATEL-ISM的采空区稳定性评价模型；依托具体的矿山采空区开展改进后的模糊综合评价法工程验证，进行采空区的稳定性定量分级。该方法能够全面考虑矿山采空区稳定性影响因素，降低传统的模糊综合评价法中的权重确定环节的主观性。本文研究成果为采空区稳定性评价提供了一种新思路，能够为金属矿山采空区的安全管理工作提供有益参考。

1 改进模糊综合评价模型

1.1 影响因素确定

通过现场调研、文献检索、问卷调查和专家咨询等方式，从岩体结构、采空区形态和环境影响等3个维度中筛选出15个影响采空区稳定性的风险因素。其中，问卷调查对象包括一线工作者、矿业工程人员（包括部门经理、车间主任和技术员等）以及金属矿山地下工程研究专家和学者，共收回100份调查问卷。通过问卷统计，得出每个风险因素的平均分数，然后再咨询行业专家，合并一些交叉的风险因素，最终确定了采空区稳定性问题涉及的3个维度共12个主要危险影响因素，如图1所示。

1.2 DEMATEL-ISM模型构建

（1）确定系统要素（F₁，F₂，…，F₁₂），构建直接影响矩阵 O，表示为

O =［x_ij ］ _m×n （1）

式中：因素x_ij （i=1，2，…，m；j=1，2，…，n；i≠j）表示因素x_i 对x_j 的直接影响；若i=j，则x_ij =0。

（2）采用行和最大值的归一化方法求得规范化矩阵 N，表示为

\mathit{N}=\mathit{O}\frac{\mathrm{1}}{\underset{\mathrm{1}\le i\le n}{\mathrm{m}\mathrm{a}\mathrm{x}} {\displaystyle \sum _{j=\mathrm{1}}^n} x_{ij}}（2）

（3）计算综合影响矩阵 T，表示为

T =［t_ij ］ _n×n （3）

综合影响矩阵 T 与规范化矩阵 N 之间的关系，表示为

\mathit{T}=\mathit{N}{(\mathit{I}-\mathit{N})}^{-\mathrm{1}}（4）

式中： I 为单位矩阵。

（4）计算各因素的影响度D、被影响度C、中心度M和原因度R。影响度是指综合影响矩阵 T 中各行值之和（Liu et al.，2023），表示该因素对所有其他因素的影响程度，记为D_i，被影响度是指综合影响矩阵 T 中各列值之和，表示该因素受其他要素影响的程度，记为C_i；中心度表示该因素在评价体系中的位置及其所起作用的大小，由要素i的影响度和被影响度相加，记为M_i，原因度表示该因素对系统内其他因素的影响程度与其受到其他因素的影响程度之间的差异，由要素i的影响度和被影响度相减得到，记为R_i。计算公式为

D_i={\displaystyle \sum _{j=\mathrm{1}}^n} x_{ij}\left(i=\mathrm{1,2},\cdots ,n\right)（5）

C_{\mathrm{i}}={\displaystyle \sum _{j=\mathrm{1}}^n} x_{ij}\left(i=\mathrm{1,2},\cdots ,n\right)（6）

M_i=D_i+C_i（7）

R_i=D_i-C_i（8）

（5）建立邻接矩阵关系矩阵 A。邻接矩阵用于表示各因素之间的直接关系，其中：a_ij =1，表示存在从x_i 到x_j 的关系；a_ij =0，表示不存在直接从x_i到x_j的关系。截距用于区分可达矩阵中的不同层级，通过设定一个阈值来消除一些间接影响，简化模型结构。

（6）建立可达矩阵 R。可达矩阵用于描述一个因素通过一系列的直接影响到达另一个因素（张玮玮等，2024），代表要素之间的传递是否会带来间接的影响关系。若数字为1，则表示某要素到另一要素之间存在路径；若数字为0，则表示某要素到另一要素之间不存在路径。可达矩阵计算过程为

{(\mathit{A}+\mathit{I})}^{k-\mathrm{1}}\ne {(\mathit{A}+\mathit{I})}^k={(\mathit{A}+\mathit{I})}^{k+\mathrm{1}}=\mathit{R}（9）

式中： I 为与 A 同阶的单位矩阵，上式中的乘法运算均为布尔代数运算。

（7）对可达矩阵进行层级划分，得出采空区稳定性影响因素的层次关系，构建对应的采空区稳定性影响因素有向层次结构图。

1.3 模糊综合评价

以影响因素层次结构图为基础，设同层次影响因素权重相同，结合专家意见确定不同层级影响因素的总权重。设权重集为A=（a₁，a₂，…，a_n ），其中a_i 为第i个因素的权重，且满足∑a_i =1。将权重集A应用于模糊综合评价权重确定环节，对原有权重确定方法进行了改进。评价等级V=｛v₁，v₂，v₃，v₄，v₅｝，分别代表很好（v₁）、良好（v₂）、一般（v₃）、较差（v₄）和非常差（v₅）。结合该采空区现场调研数据及专家分析，得到评价矩阵 W。

将权重向量 A 与评价矩阵 W 相乘，得到综合评价向量 B。在加权平均法中，每个评价等级的得分是权重与隶属度的乘积之和。根据最大隶属度原则，选择综合评价向量 B 中最大值对应的评价等级，即为采空区稳定性等级。

2 基于改进模糊综合评价的采空区稳定性评价

2.1 采空区稳定性影响因素权重分析及确定

某铅锌矿采用充填采矿法开采。其中，对于厚度小于6 m的矿体，采用平底结构浅孔留矿嗣后充填法进行回采；对于厚度为6~24 m的矿体，采用沿走向布置的分段凿岩阶段空场嗣后充填法进行回采；对于厚度大于24 m的厚大矿体，采用垂直走向布置的两步骤胶结充填采矿法。现有采空区是采用沿走向布置的分段凿岩阶段空场嗣后充填法和浅孔留矿法进行回采形成的，具体分布位置为930 m中段69线~65线之间。由于充填站建设滞后，采空区尚未充填，留设顶柱和间柱保持采空区的稳定。开展采空区稳定性评价工作，对于空区治理方案选择以及空区充填顺序确定具有促进作用。本研究选择930 m中段矿体36-2相关采场采空区进行稳定性评价，36-2采空区最大长度为90 m，最大宽度为5 m，最大高度为24 m，模型三维视图如图2所示。采空区稳定性评价技术路线如图3所示。

采用德尔菲法，邀请25名金属矿山采空区相关领域的专家学者和25名矿山企业管理工作者组成专家小组，对12个影响因素的重要程度进行两两比较。各因素相互影响程度打分情况如下：无影响（0分）、较小影响（1分）、弱影响（2分）、一般影响（3分）、较强影响（4分）和强影响（5分）。对各因素进行赋值，确定各因素间的相互影响程度，评分表形式如表1所示。通过专家评分确定不同因素之间的直接影响程度，求得的算术平均值即为每项指标的评分结果，由此构建直接影响矩阵 O，表示为

\mathit{O}=\left[\begin{array}{cccccccccccc}\mathrm{0}& \mathrm{4}& \mathrm{4}& \mathrm{3}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{2}& \mathrm{3}& \mathrm{2}\\ \mathrm{2}& \mathrm{0}& \mathrm{3}& \mathrm{2}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{2}& \mathrm{0}& \mathrm{2}& \mathrm{4}& \mathrm{2}\\ \mathrm{3}& \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{2}& \mathrm{2}& \mathrm{3}& \mathrm{3}& \mathrm{2}\\ \mathrm{3}& \mathrm{1}& \mathrm{2}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{3}& \mathrm{2}& \mathrm{0}& \mathrm{2}& \mathrm{4}\\ \mathrm{0}& \mathrm{0}& \mathrm{4}& \mathrm{1}& \mathrm{0}& \mathrm{3}& \mathrm{4}& \mathrm{4}& \mathrm{4}& \mathrm{1}& \mathrm{3}& \mathrm{3}\\ \mathrm{0}& \mathrm{0}& \mathrm{4}& \mathrm{2}& \mathrm{2}& \mathrm{0}& \mathrm{5}& \mathrm{2}& \mathrm{2}& \mathrm{3}& \mathrm{3}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{3}& \mathrm{1}& \mathrm{4}& \mathrm{3}& \mathrm{0}& \mathrm{5}& \mathrm{3}& \mathrm{2}& \mathrm{4}& \mathrm{2}\\ \mathrm{0}& \mathrm{0}& \mathrm{3}& \mathrm{3}& \mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{4}& \mathrm{4}& \mathrm{4}& \mathrm{2}\\ \mathrm{0}& \mathrm{0}& \mathrm{4}& \mathrm{2}& \mathrm{0}& \mathrm{2}& \mathrm{4}& \mathrm{2}& \mathrm{0}& \mathrm{4}& \mathrm{4}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{5}& \mathrm{2}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{4}& \mathrm{0}& \mathrm{3}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{4}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{2}& \mathrm{3}& \mathrm{3}& \mathrm{3}& \mathrm{0}& \mathrm{1}\\ \mathrm{3}& \mathrm{0}& \mathrm{4}& \mathrm{3}& \mathrm{3}& \mathrm{2}& \mathrm{2}& \mathrm{4}& \mathrm{2}& \mathrm{0}& \mathrm{3}& \mathrm{0}\end{array}\right]（10）

式中：因素x_ij （i=1，2，…，m；j=1，2，…，n；i≠j）表示因素x_i 对x_j 的直接影响；若i =j，则x_ij =0。

综合影响矩阵 T 可表示为

\mathit{T}=\left[\begin{array}{cccccccccccc}\mathrm{0.119}& \mathrm{0.192}& \mathrm{0.438}& \mathrm{0.270}& \mathrm{0.057}& \mathrm{0.062}& \mathrm{0.116}& \mathrm{0.246}& \mathrm{0.215}& \mathrm{0.274}& \mathrm{0.377}& \mathrm{0.237}\\ \mathrm{0.170}& \mathrm{0.048}& \mathrm{0.386}& \mathrm{0.224}& \mathrm{0.055}& \mathrm{0.060}& \mathrm{0.113}& \mathrm{0.262}& \mathrm{0.207}& \mathrm{0.263}& \mathrm{0.386}& \mathrm{0.220}\\ \mathrm{0.197}& \mathrm{0.084}& \mathrm{0.297}& \mathrm{0.197}& \mathrm{0.057}& \mathrm{0.066}& \mathrm{0.122}& \mathrm{0.261}& \mathrm{0.278}& \mathrm{0.303}& \mathrm{0.362}& \mathrm{0.218}\\ \mathrm{0.218}& \mathrm{0.091}& \mathrm{0.402}& \mathrm{0.188}& \mathrm{0.081}& \mathrm{0.087}& \mathrm{0.181}& \mathrm{0.339}& \mathrm{0.305}& \mathrm{0.230}& \mathrm{0.372}& \mathrm{0.313}\\ \mathrm{0.152}& \mathrm{0.058}& \mathrm{0.652}& \mathrm{0.309}& \mathrm{0.137}& \mathrm{0.256}& \mathrm{0.394}& \mathrm{0.511}& \mathrm{0.529}& \mathrm{0.406}& \mathrm{0.574}& \mathrm{0.368}\\ \mathrm{0.132}& \mathrm{0.052}& \mathrm{0.581}& \mathrm{0.298}& \mathrm{0.185}& \mathrm{0.128}& \mathrm{0.381}& \mathrm{0.393}& \mathrm{0.414}& \mathrm{0.415}& \mathrm{0.507}& \mathrm{0.270}\\ \mathrm{0.144}& \mathrm{0.055}& \mathrm{0.617}& \mathrm{0.305}& \mathrm{0.262}& \mathrm{0.250}& \mathrm{0.259}& \mathrm{0.532}& \mathrm{0.500}& \mathrm{0.433}& \mathrm{0.597}& \mathrm{0.333}\\ \mathrm{0.126}& \mathrm{0.049}& \mathrm{0.511}& \mathrm{0.315}& \mathrm{0.122}& \mathrm{0.139}& \mathrm{0.226}& \mathrm{0.280}& \mathrm{0.442}& \mathrm{0.417}& \mathrm{0.496}& \mathrm{0.278}\\ \mathrm{0.123}& \mathrm{0.049}& \mathrm{0.545}& \mathrm{0.280}& \mathrm{0.101}& \mathrm{0.176}& \mathrm{0.321}& \mathrm{0.357}& \mathrm{0.314}& \mathrm{0.424}& \mathrm{0.502}& \mathrm{0.245}\\ \mathrm{0.091}& \mathrm{0.037}& \mathrm{0.422}& \mathrm{0.201}& \mathrm{0.047}& \mathrm{0.061}& \mathrm{0.117}& \mathrm{0.178}& \mathrm{0.322}& \mathrm{0.184}& \mathrm{0.331}& \mathrm{0.168}\\ \mathrm{0.095}& \mathrm{0.038}& \mathrm{0.440}& \mathrm{0.198}& \mathrm{0.069}& \mathrm{0.081}& \mathrm{0.203}& \mathrm{0.310}& \mathrm{0.336}& \mathrm{0.323}& \mathrm{0.279}& \mathrm{0.192}\\ \mathrm{0.253}& \mathrm{0.074}& \mathrm{0.612}& \mathrm{0.360}& \mathrm{0.216}& \mathrm{0.200}& \mathrm{0.296}& \mathrm{0.479}& \mathrm{0.426}& \mathrm{0.338}& \mathrm{0.535}& \mathrm{0.260}\end{array}\right]（11）

根据式（4）~式（8）计算各度量值，计算结果如表2所示。

对中心度归一化处理即得到各影响因素的权重，分别为：F₁（5.78%）、F₂（4.21%）、F₃（10.91%）、F₄（7.78%）、F₅（7.50%）、F₆（6.96%）、F₇（9.17%）、F₈（9.87%）、F₉（10.10%）、 F₁₀（8.06%）、F₁₁（10.30%）和F₁₂（9.35%）。

综合影响矩阵的平均值为0.266，标准差为0.152，确定截距λ为两项之和0.418，若t _ij ≥λ，则t _ij =1；若t_ij ≤λ，则t_ij =0，得到布尔邻接矩阵 A。由式（9）求得可达矩阵 R 为

\mathit{R}=\left[\begin{array}{cccccccccccc}\mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{1}& \mathrm{1}\end{array}\right]（12）

按照结果优先的方式对可达矩阵进行层级划分，构建出对应的采空区稳定性影响因素有向层次结构图（图4），由图可知采空区稳定性影响因素的层次关系与重要程度。

该结构图由10个带边的节点、2个孤立的节点和19条有向边组成，代表着不同风险因素之间的影响关系。每个节点代表一个采空区稳定性事故的影响因素。由图4可知，采空区稳定性影响因素具有多级递进的特点，由下至上划分为5个层次。节点之间的有向边表示2个风险因素之间的影响关系。结合DEMATEL-ISM模型，对采空区事故的复杂关系进行分析，具体如下：

（1）根据DEMATEL方法，表2显示了采空区事故中各危险因素的影响度、被影响度、原因度和中心度。其中，原因度较高的因素，如空区体积（F₅）和空区跨度（F₇）等对其他因素产生影响，因此判定为主要影响因素。中心度较高的因素，如空区时间效应（F₈）和支护条件（F₁₁）等容易受其他危险因素的影响。对于管理层来说，重要的是要防止其他可能导致系统性风险因素的干扰。在实际操作中，建议从本质安全的角度优先关注原因因素。如：空区体积（F₅）被确定为原因因素，该因素对空区稳定性具有直接的重大影响，空区体积参数的变化可引发连锁反应，导致一系列事故。

（2）根据ISM方法，图5中将采空区事故的危险因素系统划分为本质原因、过渡原因和直接原因3个层次。第一级别L₁代表直接风险，这是造成事故的最直接因素。L₂至L₄级代表过渡风险，起着桥梁作用。L₅级代表本质风险，能够从根本上降低事故发生的风险。该结论为安全管理人员提供了控制有效性和安全投资方面的建议。

2.2 采空区稳定性评价结果

确定综合风险重要性计算公式为

S=L_{\mathrm{5}}\times \mathrm{40}\mathrm{\%}+L_{\mathrm{4}}\times \mathrm{30}\mathrm{\%}+L_{\mathrm{3}}\times \mathrm{20}\mathrm{\%}+L_{\mathrm{2}}\times \mathrm{20}\mathrm{\%}+L_{\mathrm{1}}\times \mathrm{10}\mathrm{\%}（13）

式中：S为综合风险等级；L_i 为采空区稳定性影响因素层次结构图中第i层的影响因素的重要程度。将每层中的影响因素权重视为相等，则权重集A=（0.040，0.033，0.033，0.033，0.033，0.133，0.133，0.300，0.040，0.040，0.040，0.133）。

结合该采空区现场调研数据及专家分析，得到以下评价矩阵 W，其中w_ij 表示因素a_i 对评价等级v_j 的隶属度。即所评价空区的12个稳定性影响因素对评价等级的隶属程度。

\mathit{W}=\left[\begin{array}{ccccc}\mathrm{0.1}& \mathrm{0.3}& \mathrm{0.3}& \mathrm{0.2}& \mathrm{0.1}\\ \mathrm{0.2}& \mathrm{0.2}& \mathrm{0.2}& \mathrm{0.2}& \mathrm{0.2}\\ \mathrm{0.1}& \mathrm{0.2}& \mathrm{0.3}& \mathrm{0.2}& \mathrm{0.1}\\ \mathrm{0.3}& \mathrm{0.4}& \mathrm{0.1}& \mathrm{0.1}& \mathrm{0.1}\\ \mathrm{0.1}& \mathrm{0.2}& \mathrm{0.3}& \mathrm{0.2}& \mathrm{0.2}\\ \mathrm{0.3}& \mathrm{0.2}& \mathrm{0.2}& \mathrm{0.2}& \mathrm{0.1}\\ \mathrm{0.1}& \mathrm{0.1}& \mathrm{0.2}& \mathrm{0.3}& \mathrm{0.3}\\ \mathrm{0.2}& \mathrm{0.3}& \mathrm{0.2}& \mathrm{0.2}& \mathrm{0.1}\\ \mathrm{0.2}& \mathrm{0.3}& \mathrm{0.4}& \mathrm{0.1}& \mathrm{0.0}\\ \mathrm{0.3}& \mathrm{0.3}& \mathrm{0.2}& \mathrm{0.1}& \mathrm{0.1}\\ \mathrm{0.4}& \mathrm{0.4}& \mathrm{0.1}& \mathrm{0.1}& \mathrm{0.0}\\ \mathrm{0.2}& \mathrm{0.4}& \mathrm{0.1}& \mathrm{0.1}& \mathrm{0.0}\end{array}\right]（14）

然后，求出综合评价向量 B，表示为

\mathit{B}=\mathit{A}\times \mathit{W}=\{\mathrm{0.203,0.268,0.196,0.183,0.111}\}（15）

因此，综合评价向量 B =（B₁，B₂，B₃，B₄，B₅）=（0.203，0.268，0.196，0.183，0.111）。 B₂的值最大，为0.268，因此确定评价等级为v₂（良好）。由此可知，930 m中段36-2采空区的稳定性等级为良好。

3 评价结果检验

在完成采空区稳定性评价之后，采用数值模拟方法（Li et al.，2024）进行采空区稳定性数值模拟验证，从而评价结果的准确性。

3.1 FLAC^3D数值分析模型建立

利用空区三维扫描获得的几何数据以及矿山相关资料，针对930 m中段、870 m中段、820 m中段和历史采空区，构建了全矿尺度的FLAC^3D计算模型，如图5所示。模型尺寸为2 300 m×2 050 m×（680~968）m，该模型共包含约650万个单元，并考虑了现场扫描的巷道和空区开挖空间，基于测绘资料对重构的采空区模型进行修正。

在模拟过程中，岩体采用应变软化本构模型（沈华章等，2016），充填体采用摩尔—库伦本构模型，开挖采用空模型进行计算。利用SketchUp软件提取矿区三维地表地形图和地表起伏数据，建立该矿区的地表三维模型。

3.2 模拟参数确定

通过现场点荷载以及室内岩石力学单轴和三轴试验，得到矿区3个中段上盘、矿体和下盘岩石的变形和强度参数，采用Hoek-Brown准则估算对应的岩体物理力学参数（张景轩，2022）。通过计算，得到上盘、矿体和下盘岩体的变形和强度参数如表3所示。在后续的模拟中，采用类似试验数据得出的强度和弹模关系作为经验公式。

在模拟过程中，充填体的泊松比取值为0.30，内摩擦角取值为33°。本模型地应力分布情况如式（16）所示（何建，2017）。此次模拟的3个中段主要集中在14~250 m区间，属于典型的浅部埋深情况（黄庆享，2002）。

\left\{\begin{array}{l}{\sigma }_V=\mathrm{0.0271}\times H\\ {\sigma }_{\mathrm{1}}=\mathrm{0.0216}\times H+\mathrm{6.7808}\\ {\sigma }_{\mathrm{2}}=\mathrm{0.0182}\times H+\mathrm{2.2328}\end{array}\right.（16）

3.3 采空区顶柱稳定性分析结果

图6显示了相关采场采空区顶部岩体差应力分布状态。由于采空区顶板距离山体约为152 m，且近似平行于最大水平主应力方向，岩体发生的破坏类型为拉破坏，故在采空区周边有一定的压力升高区，局部可达20 MPa（{\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}），在周边围岩体内也有一定范围的拉应力屈服和少量的剪切破坏发生。采空区上部围岩的屈服极限差应力最大值为15 MPa（{\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}），根据第三强度理论，有以下关系式：

\frac{\mathrm{1}}{\mathrm{2}}{\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}\le ({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}})（17）

最大剪应力不超过屈服极限，未达到破坏条件，表明36号矿体的上覆岩层较稳定，不易发生大范围顶板冒落。当前，采空区处于稳定状态，相应的采空区顶板与周围岩体也处于稳定状态。因此，数值分析结果与本文提出的采空区稳定性评价结果（良好）一致，验证了采空区稳定性评价模型的准确性。

3.4 现场空区概况

根据矿山现场观察，930 m中段36-2采空区地质结构坚实，不存在大的断裂带或软弱层，岩石力学性质稳定，采空区周围的围岩没有明显的变形或破坏迹象，能够维持原有的力学性能；开采计划也充分考虑了地质条件，采取合理的采动顺序和开采方法，确保了采空区的稳定性，该空区不会发生大规模的塌陷或位移。间隔半年后对空区进行三维激光扫描形态对比（图7），空区形态无明显变化，与评价结果相吻合。

4 结论

（1）通过专家调查确定了12个与采空区事故相关的风险因素，由此组成风险因素集，该因素集尽可能多地涵盖了采空区稳定性影响因素。采用DEMATEL-ISM模型作为风险分析方法，确定了事故影响因素的因果属性和层次关系，并根据分析结果给出了各影响因素的权重，改进了模糊综合评价中影响因素权重获取方式。

（2）提出了一种改进权重确定的模糊评价方法用于采空区稳定性风险评价，该方法能够实现采空区稳定性风险评价及分级，为采空区稳定性评价提供了理论依据和技术支撑。将所提方法应用于某铅锌矿930 m中段36-2采空区，确定该空区风险等级为二级，评价结果与数值模拟和现场测试结果相吻合，表明所提方法能够系统、全面地反映各采空区稳定性影响因素的风险程度。

（3）采用FLAC^3D数值模拟软件，从相关采场采空区顶部岩体的差应力来模拟判断采空区顶板稳定性状况，分析结果表明该采空区上覆岩层比较稳定，采空区稳定性良好，与模糊评价方法对采空区稳定性分级结果基本一致，验证了所提评价方法对采空区稳定性分级的准确性。

参考文献

[1] Cai Xin， Huang Yashuai， Zhou Zilong，et al，2024.Redefinition and identification of critical pillars in goaf based on the risk analysis of pillar instability［J］.Chinese Journal of Rock Mechanics and Engineering，43（1）：216-225.

[2] Cai Yong， Deng Guoliang， Zhen Libing，et al，2022.Mechanical model analysis of goaf considering pillar rheological effect［J］.Metal Mine，51（2）：96-101.

[3] Dalvi-Esfahani M， Niknafs A， Kuss D J，et al，2019.Social media addiction：Applying the DEMATEL approach［J］.Telematics and Informatics，43（Oct.）：101250.

[4] Guo Huimin， Cheng Lianhua， Li Shugang，2023.Research on causal factors of coal mine gas explosion based on DEMATEL-ISM-MICMAC［J］.Mining Safety and Environmental Protection，50（2）：114-119.

[5] Guo Q， Li Y， Meng X，et al，2019.Instability risk assessment of expressway construction site above an abandoned goaf：A case study in China［J］.Environmental Earth Sciences，78：1-14.

[6] He Biaoqing，2018.Numerical simulation on stability of large goaf group base on FLAC^3D ［J］.China Mining Magazine，27（7）：99-102.

[7] He Jian，2017.The Study of Geostress Characteristics and Engineering Areas Geostress Inversion in the Southwest Region［D］.Chongqing：Chongqing University.

[8] Hong Weibin， Sheng Wu，2022.Analysis of influencing factors of coal mine water inrush accidents based on DEMATEL-ISM-BN［J］.Journal of Mine Automation，48（12）：116-122.

[9] Hu Y X， Li X B，2012.Bayes discriminant analysis method to identify risky of complicated goaf in mines and its application［J］.Transactions of Nonferrous Metals Society of China，22（2）：425-431.

[10] Huang Qingheng，2002.Ground pressure behavior and definition of shallow seams［J］.Chinese Journal of Rock Mechanics and Engineering，21（8）：1174-1177.

[11] Jiang L C， Jiao H Z， Wang Y D，et al，2021.Comprehensive safety factor of roof in goaf under deep high stress［J］.Journal of Central South University，28（2）：595-603.

[12] Ke Lihua， Meng Yaoyao， Yao Nan，et al，2023.Evaluation of goaf stability based on fuzzy statistical method［J］.China Safety Science Journal，33（2）：59-67.

[13] Li Junping， Zhou Chuangbing， Li Xiangyang，2005.Similarity simulation of disposal scheme for abandoned stope under the sunken topography［J］.Chinese Journal of Rock Mechanics and Engineering，（4）：581-586.

[14] Li M L， Li K G， Liu Y D，et al，2024.Goaf risk prediction based on IAOA-SVM and numerical simulation：A case study［J］.Underground Space，15（2）：153-75.

[15] Liu T， Ding C， Feng H，et al，2023.Risk evaluation for oil tanker berthing and handling operations using the improved DEMATEL-ISM approach［J］.Journal of Loss Prevention in the Process Industries，86：105-182.

[16] Shen Huazhang， Wang Shuilin， Guo Mingwei，et al，2016.A preliminary study of the progressive failure and stability of slope with strain-softening behaviour［J］.Rock and Soil Mechanics，37（1）：175-184.

[17] Wang F， Zou P， Meng Z，et al，2019.Study on stability of goaf pillars in Daqiao phosphate mine：Theoretical calculation and field investigation［J］.Geotechnical and Geological Engineering，37：1483-1492.

[18] Wang L， Guo Q， Yu X，2023.Stability-level evaluation of the construction site above the goaf based on combination wei-ghting and cloud model［J］.Sustainability，15（9）：7222.

[19] Warfield J N，1974.Developing interconnection matrices in structural modeling［J］.IEEE Transactions on Systems，Man and Cybernetics，（1）：81-87.

[20] Xie Lianku， Xiong Daiyu， Yang Tianhong，et al，2020.Analysis on disaster impact of concealed goaves underneath open-pit slope and its treatment［J］.The Chinese Journal of Nonferrous Metals，30（10）：2503-2512.

[21] Zhang Jingxuan，2022.Study on Overburden Structure Characteristicsand Gob-side Entry Retention Technology of Lower Stratied Working Face in Guhanshan Mine［D］.Jiaozuo：Henan Polytechnic University.

[22] Zhang Q L， Wang C H， Han L，et al，2023.Study on stability ana-lysis and treatment of underground goaf in metal mine［J］.Mining，Metallurgy and Exploration，40（5）：1973-1985.

[23] Zhang Weiwei， Yang Chun， Yao Kaiwen，et al，2024.Research on livelihood strategy selection of resettlement based on DEMATEL-ISM model［J］.China Rural Water and Hydropower，（6）：162-169.

[24] Zhang X， Hu B， Zou J，et al，2021.Quantitative characterization of overburden rock development pattern in the goaf at different key stratum locations based on DEM［J］.Advances in Civil Engineering，2021：8011350.

[25] Zhao B， Zhao Y， Wang J，2021.New stability forecasting model for goaf slope based on the AHP-TOPSIS theory［J］.Arabian Journal of Geosciences，14：1-12.

[26] Zhou Keping， Cao Lixiong， Li Jielin，et al，2022.Stability numerical analysis and safety classification evaluation of complex goaf group［J］.Gold Science and Technology，30（3）：324-332.

[27] 蔡鑫，黄涯帅，周子龙，等，2024.基于矿柱失稳风险分析的采空区关键矿柱重定义与辨识［J］.岩石力学与工程学报，43（1）：216-225.

[28] 蔡泳，邓国梁，甄利兵，等，2022.考虑矿柱流变效应的采空区力学模型分析［J］，金属矿山，51（2）：96-101．

[29] 郭慧敏，成连华，李树刚，2023.基于DEMATEL-ISM-MICMAC的煤矿瓦斯爆炸致因研究［J］.矿业安全与环保，50（2）：114-119．

[30] 何标庆，2018.基于FLAC^3D的大型采空区群稳定性数值模拟分析［J］.中国矿业，27（7）：99-102．

[31] 何建，2017.西南地区地应力特征及工程区域地应力反演研究［D］.重庆：重庆大学.

[32] 洪伟斌，盛武，2022.基于DEMATEL-ISM-BN的煤矿透水事故影响因素分析［J］.工矿自动化，48（12）：116-122．

[33] 黄庆享，2002.浅埋煤层的矿压特征与浅埋煤层定义［J］.岩石力学与工程学报，21（8）：1174-1177.

[34] 柯丽华，孟瑶尧，姚囝，等，2023.基于模糊统计法的采空区稳定性评价［J］.中国安全科学学报，33（2）：59-67．

[35] 李俊平，周创兵，李向阳，2005.下凹地形下采空区处理方案的相似模拟研究［J］.岩石力学与工程学报，（4）：581-586．

[36] 沈华章，王水林，郭明伟，等，2016.应变软化边坡渐进破坏及其稳定性初步研究［J］.岩土力学，37（1）：175-184．

[37] 解联库，熊代余，杨天鸿，等，2020.露天采场边坡下隐伏空区灾害影响分析及其治理［J］.中国有色金属学报，30（10）：2503-2512.

[38] 张景轩，2022.古汉山矿下分层工作面覆岩结构特征及沿空留巷技术研究［D］.焦作：河南理工大学.

[39] 张玮玮，杨春，姚凯文，等，2024.基于DEMATEL-ISM模型的水库移民生计策略选择研究［J］.中国农村水利水电，（6）：162-169．

[40] 周科平，曹立雄，李杰林，等，2022.复杂采空区群稳定性数值分析及安全分级评价［J］.黄金科学技术，30（3）：324-332.

基金资助

国家自然科学基金项目“热—应力—腐蚀耦合作用下自膨胀锚杆锚固强度时效性能研究”(52074062)

“爆破载荷下自膨胀填塞体作用机理与爆破能量调控机制研究”(51874068)