地下金属矿山采用二步骤开采时会留设大量起支撑作用的矿柱，形成矿柱群。随着开采的不断推进，岩体原有应力平衡被打破，矿柱群所承受的应力将是原岩应力的数倍。当矿柱群所受应力超过其强度时就会发生失稳破坏，引发顶板冒落、矿柱垮塌、岩层移动、地表崩落和大范围垮塌冲击等现象（李锋，2014；何荣兴等，2020）。矿柱群失稳不仅会扰乱正常的生产秩序，而且会严重威胁矿区作业人员的生命安全（赵康等，2016；陈顺满等，2018；叶光祥等，2021）。因此，如何对矿柱群进行稳定性分析，计算出矿柱群失稳概率，采取科学有效的防治措施进行矿柱群治理并安全回采矿柱是当前亟待解决的重要技术难题。

针对该问题，国内外学者采用数值模拟、位移监测和概率分析等方法进行了大量研究。熊有为等（2023）提出了隐伏空区群综合探测及数字化建模技术，用来分析隔离矿柱稳定性；周子龙等（2018a）构建了并联矿柱承载模型，对矿柱在不同安全系数下的稳定性进行了评价；Cording et al.（2015）通过数值模拟研究了单个矿柱失稳引起整个矿柱群失稳的动态过程。此外，鉴于矿柱群失稳灾变机理的复杂性，Najafi et al.（2011）采用概率法克服了矿柱群失稳研究中不确定性带来的影响。以上研究表明，针对单个矿柱的稳定性分析已取得许多重要成果，但多集中于因本身承载能力不足导致单个矿柱的失稳过程，而对多矿柱间荷载传递导致整体矿柱群失稳的研究比较欠缺。

矿柱群失稳是一个自组织过程（罗一忠，2005；谢杰辉，2018），通常情况下都是单个矿柱内小裂隙的拓展贯通，随着一系列的发展演化成大裂隙直至矿柱破裂失稳，单个矿柱破裂后荷载传递到相邻矿柱又引起相邻矿柱失稳，随着失稳的单个矿柱数量增多，最终导致整体矿柱群级联失稳。周子龙等（2018b）基于PFC模拟了矿柱大规模级联失稳的发生机理及载荷传递过程；丘帆等（2014）基于Voronoi图预测了不同时间段的矿柱群失稳情况；Ma et al.（2012）研究了矿柱发生大规模级联失稳的“多米诺效应”。从系统学观点考虑，整个矿柱群失稳过程存在一个临界状态，因此可基于重整化群理论对该过程给予准确的描述。李恩翀等（2021）基于一维重整化群模型计算出矿柱失稳概率，并通过与数值模拟结果进行比对验证了该模型的准确性；朱德福等（2018）计算了采空区中煤柱群一维、二维重整化群模型对应的临界失稳概率；张淑坤等（2016）研究了重整化群模型中不同垂直应力系数下煤柱群临界失稳概率的大小。以上研究表明，针对单个矿柱的失稳以及多矿柱级联失稳已取得许多重要的成果。然而，这些研究多集中于对矿柱群失稳机理及发生的动态过程的研究，而缺少对单个矿柱失稳导致矿柱群整体失稳的定量评估研究。

因此，本文基于数学方法对单个矿柱的失稳风险进行了定量评估，研究了单个矿柱失稳导致矿柱群失稳的概率变化，为下一步确定最低风险的矿柱群回采路线提供理论指导。以破碎盘区矿柱群为研究对象，从重整化群理论（RG）和贝叶斯网络（BN）2个方面建立稳定性分析模型，对铜坑矿矿柱群稳定性进行分析，通过概率更新为矿柱安全回采方案的选择提供决策依据，并为同类型矿山矿柱群稳定性评价及安全回采提供经验借鉴。

矿柱群失稳是一个渐进破坏的过程。如图1（a）所示，未开采前岩体处于稳定状态。图1（b）表示矿柱群形成初期，随着应力的重新分布，高应力集中储存在矿柱中，此时矿柱受力处于弹性阶段，根据弹塑性理论，当矿柱受力处于弹性阶段时，荷载移除后的弹性变形可以恢复，矿柱群系统会保持长期稳定状态。图1（c）表示随着时间的推移，受采场开采活动、顶板荷载增加和矿柱强度降低等因素的影响，矿柱强度达到峰值。图1（d）表示矿柱中的高应变能开始释放，局部矿柱向失稳趋势发展；在损伤进一步发展至核心承载区时，矿柱的支撑作用消失，原有的应力平衡被打破，该区域的顶板会发生局部冒落形成新的应力拱，为适应新的平衡状态，原来由采空区位置岩体承担的荷载转移至周边矿柱及围岩体（张海波等，2014），局部矿柱开始失稳。图1（e）表示随着矿柱应力的传递，应力相对集中部位的应变能通常会转移至相邻强度最低的单元矿柱，周边矿柱承受的荷载增加，该区域应力集中程度增高。当所承受的应力超过其临界值时，又会导致临近矿柱的进一步失稳，矿柱中的弹性能大量释放，该区域顶板应力又将传递给次一级强度较低的单元矿柱，形成如图1（f）所示的类似“多米诺骨牌”效应（朱万成等，2022），最终当足够多的矿柱发生失稳时，矿柱群整体失稳，此时采空区将发生大面积瞬时塌陷。

矿柱失稳多发生于高压应力集中的区域，开采深度的增加导致压应力不断增加，采场中的裂纹不断扩展延伸，裂纹扩展的速率和范围逐渐由稳态向非稳态转变。当矿柱的稳定性处于临界状态时，任何微小的扰动（微扰）都不可忽视（王中秋，2012）。现代非线性科学认为，在临界点处，扰动的诱发作用主要是通过扰动的放大效应来实现（张燕等，2013）。由于这时系统处于高度不稳定状态，任何微小的扰动都会被放大，微扰动在临界点附近会转变成巨扰动，正是这种巨扰动驱动着事物向新的状态演化，所以在一定情况下裂纹非稳定扩展的起始点可视为应变能快速释放的开端，最终可视为矿柱失稳发生的临界点。这种现象有一定的自组织性，重整化群理论（RG）可以很好地解释该现象，重整化群能从小的单元合并成大的单元，在这一过程中，系统小范围的涨落被逐渐消除，从而突出与长程关联的因素，整个过程并不会改变系统的物理性质。

贝叶斯网络（BN）是一种基于概率图形的网络，又称为信度网络，它是一种有向无环图模型，在处理信息不完全和不确定性的小样本事件时具有明显优势，且随着信息的更新可以实现随着时间发展的动态概率评估（张晓全等，2009）。近年来，该方法在矿山灾害治理领域得到了广泛应用（胡建华等，2022；李敏等，2023）。

本研究基于贝叶斯网络构建了矿柱群失稳概率模型（图2），模型主要由节点、有向弧和条件概率表3个部分组成。子节点代表单个矿柱失稳事件，父节点代表矿柱群整体失稳事件，节点之间用有向弧连接。弧表示单个矿柱失稳与矿柱群整体失稳的逻辑关系，条件概率表表示连接关系强度，其中S、U和F分别代表矿柱稳定、不稳定和失稳概率。基于贝叶斯网络构建的矿柱群失稳概率模型，能够充分利用单个矿柱形态和岩石力学参数，更直观地表达单个矿柱失稳与矿柱群整体失稳的概率关系。图形化的知识表示方法使得保持概率知识库的一致性和完整性变得容易，当某些因素发生改变时，便于进行贝叶斯网络模块的动态更新（李金蓉等，2022）。

基于已有研究，对破碎盘区矿柱群稳定性展开研究。通过综合考虑重整化群理论（RG）和贝叶斯网络（BN）的优势，构建了基于RG-BN理论的破碎盘区矿柱群稳定性分析模型。矿柱群稳定性分析流程如图3所示。

基于二维RG方法进行矿柱群失稳临界概率分析。根据重整化群理论，首先建立如图4所示的二维RG模型，并假定图中黑色单元为已失稳状态，白色单元为未失稳状态。其中，蓝色区域代表一级原胞，二维RG模型下共有4个基本单元；红色区域代表二级原胞，共有16个基本单元，依此类推，直到N级原胞。图中总原胞数量记为M，已失稳原胞数量记为m，已失稳原胞占总原胞数量的比例代表模型中原胞的破坏率，即矿柱群的失稳概率。

为了便于RG模型的拓展分析，需做出如下假设：（1）每个单元的面积及所受荷载均相等，当单元所受荷载超过其单元极限承载强度临界值时，认为该单元发生失稳破坏；（2）假设基本单元破坏的概率为P₀，二维情况下一级原胞由4个单元体组成，在此研究一级原胞的失稳概率P₁。在二维情况下，一级原胞中基本单元破坏形式及概率分布情况如图4（b）所示，图中黑色部分为已失稳单元，白色为未失稳单元。

如图4（b）所示，原胞a由4个未失稳单元组成，原胞失稳概率为0；原胞b存在3个未失稳单元，1个已失稳单元，由于荷载传递，已失稳单元也可造成未失稳单元的破坏，造成该原胞失稳破坏有4种分布，失稳概率为4P0(1-P0)3Pab3（Pab为一个单元失稳后向邻近单元传递应力的条件概率）；原胞c存在2个未失稳单元，2个已失稳单元，造成该原胞失稳破坏有6种分布，失稳概率为6P02(1-P0)2Pab2；原胞d存在1个未失稳单元，3个已失稳单元，造成该原胞失稳破坏有4种分布，失稳概率为4P03(1-P0)Pab；原胞e中所有单元均已失稳，则其原胞失稳概率为P04。

对以上所有分布进行求和，则一级原胞失稳破坏概率为

考虑失稳单元对相邻未失稳单元传递荷载的影响，Pab为一个单元失稳后向邻近单元传递应力的条件概率，概率公式为

式中:Pa为单元未破坏前相邻单元的失稳概率；Pb为单元破坏后相邻单元的失稳概率。

基本单元失稳引起系统失稳的破坏概率均服从Weibull分布，其概率密度函数可表示为

可靠度函数可表示为

式中：m为关于单元强度的相关性参数；F为此时单元承载作用力；Fm为单元承载的强度极限。

将式（5）代入式（1），导出一级原胞失稳的概率表达式为

当m=3时（朱德福，2018），对式（6）进行重整化，可得N级原胞失稳的概率表达式为

令Pn=Pn-1，在0≤Pn-1≤1范围内，由Pn与Pn-1的RG映射曲线（图5），可解得3个不动点，即0、P*和1，忽略0、1这2个无意义点，则N级原胞失稳的临界概率为P*=0.3088。结合可靠度分析，当N级原胞承载力Fn=0时，代入式（3）和式（4），解得PFn=0，RFn=1，即此时N级原胞失稳概率为0，可靠度为100%；同理，当N级原胞承载力Fn等于其极限承载力Fm时，解得PFn=0.632，RFn=0.368，此时N级原胞失稳概率为0.632，可靠度为36.8%，即矿柱群极限临界失稳概率Pm=0.368,此时矿柱处于极限失稳状态。综上可知，P*可将系统划分为3个状态区域，当0≤P*≤0.3088时，经过若干次迭代，P*→0，此时N级原胞不发生失稳或仅局部发生失稳，整个系统仍然处于稳定状态；当0.3088≤P*≤0.368时，N级原胞开始向失稳趋势演化；当P*≥0.368时，经过若干次迭代，P*→1，此时N级原胞已完全失稳。

铜坑矿92号矿体为一厚大缓倾斜矿体，上部为细脉带和91号矿体，并与91号矿体部分重叠。矿体分7个盘区进行回采，盘区宽度为90~100 m，盘区矿柱沿南北向勘探线布置，盘区矿柱宽度为20 m。盘区内矿块划分为矿房和矿柱，沿矿体走向布置。矿房宽度为25 m，长度为70~80 m，高度为10~60 m；采用Φ165mm、Φ110mm大孔和Φ90mm、Φ65mm中深孔凿岩，铲运机出矿，采场生产能力为3 000~4 000 t/d。

92号矿体规模化开采近20年，由于前期主要采用崩落法回采，造成地表塌陷，应力集中使得余下的大部分盘区矿柱和房间柱较破碎，回采难度大。近年来，盘区间柱资源的回采矿石量在矿山总年产量中的占比逐渐提高，总计矿量超过600万t。由于目前矿山面临自上而下多中段同时回采的复杂局面，随着时间的推移，地压活动日益加剧，矿柱资源回采难度持续升高，亟需对该区域开展稳定性分析，并研究最安全的回采方案。

实践证明，合理的回采顺序在一定程度上能够弱化矿柱中的应力集中程度，有效释放岩体开挖聚集的能量，既能提升矿柱群的稳定性，也能确保矿柱回采的效率及安全性，预防矿柱回采过程中矿柱群大规模垮塌事件的发生。因此，92号矿体矿柱群整体回采思路为：通过稳定性分析，选取最安全回采方案，先采部分矿体时保证后续矿体的回采安全性。图6所示为该矿柱群待回采区域矿柱布局，为便于分析，按矿柱位置分布分别规划4个待回采区。其中，T209-210、T210-211、T211-212和T212-213矿柱为1号待回采区；T111-112和T115-S1矿柱为2号待回采区； T313-314、T501-502和T502-503矿柱为3号待回采区；1号、2号和3号间柱为4号待回采区。

对矿柱群进行稳定性分析的重点在于比较实际失稳概率与临界失稳概率，将单个矿柱的稳定性情况量化后输入贝叶斯网络，可计算得到矿柱群的实际失稳概率。表1为该矿柱群的矿柱形态和岩石力学参数。

单个矿柱的稳定性情况一般用安全系数表示。然而，由于岩石的各向异性、矿柱形态不规则和不同采矿方法的影响，使用确定性方法计算的安全系数具有一定的局限性，计算得到安全系数大于1的矿柱仍可能失稳，安全系数小于1的矿柱也可能稳定（黄敏等，2020）。为了更精确地识别矿柱稳定性情况，Wattimena（2014）综合考虑矿柱横断面所受的平均应力与单个矿柱极限承载强度之比和矿柱宽高比，通过多项逻辑回归模型对矿柱稳定性情况进行进一步分析，将矿柱的稳定性划分为稳定概率、不稳定概率和失稳概率3种情况。

单个矿柱横断面所受的平均应力（熊晓勃等，2021）计算公式为

式中：σP为横断面所受的平均应力（MPa）；S为矿柱承受载荷面积（m²）；g为重力加速度，取9.8m/s2；ρi为第i层上覆岩层的平均密度（kg/m³）；Hi为第i层上覆岩层的厚度（m）；l和b分别为矿柱的长度（m）和宽度（m）。

矿柱的极限承载强度计算公式为

式中：σR为矿柱的极限承载强度（MPa）；σc为岩石的单轴抗压强度。

根据矿柱稳定性数据建立多项逻辑回归模型，推导出式（10）~式（12），计算出矿柱稳定概率Ps、不稳定概率Pu和失稳概率Pf。

将表1中的各项参数代入式（10）~式（12）中，计算得到单个矿柱稳定性概率分布（表2）。

将该矿柱群单个矿柱转换为BN网络节点，单个矿柱失稳概率输入Genie软件进行参数学习，建立如图7所示的BN模型。由BN模型计算得到矿柱群失稳概率Pf=14.85%，失稳概率小于矿柱群临界失稳概率P*，判断此时矿柱群处于稳定状态。

贝叶斯网络的一个重要优势是可以进行概率反演，即在顶上事件确定发生的情况下推断出各基本因素对顶上事件的重要程度。在本模型中通过定义单个待采区失稳概率为1，即在单个待采区确定失稳的情况下，通过矿柱群失稳贝叶斯网络模型计算得到矿柱群失稳概率。图8所示为4个待采区已回采时矿柱群失稳的概率。若某个待采区回采导致矿柱群失稳概率明显增高，则表明该待采区在矿柱群失稳过程中更为敏感，反之，对矿柱群失稳影响较小的待采区即为首先选取的安全回采区域。

分析图8可知，第一步回采导致矿柱群承压矿柱减少，其余待采区发生应力重新分布和动态变化，出现了中心地段承载矿柱高应力集中现象，单个矿柱承受荷载的增加导致其失稳概率增加，进而导致矿柱群失稳概率整体增加，但回采1号、2号、3号待采区导致矿柱群失稳概率均小于矿柱群失稳临界概率P*，表明首次回采该区域并不会导致矿柱群失稳。其中，1号待采区失稳导致矿柱群失稳概率为25.51%，增幅为10.66%；2号待采区失稳导致矿柱群失稳概率为27.46%，增幅为12.61%；3号待采区失稳导致矿柱群失稳概率为24.69%，增幅为9.84%。为此，选取3号待采区作为首次回采对象，4号待采区为盘区间柱。若将间柱作为首次回采对象，则矿柱群失稳概率高达63.16%，增幅为48.31%，远超矿柱群临界失稳概率P*，必定导致矿柱群失稳，故不能将其作为首要回采对象。

图9显示，在首次回采的基础上，二次回采1号待采区导致矿柱群失稳的概率为38.59%；回采2号待采区导致矿柱群失稳概率为37.2%，回采1号和2号待采区导致矿柱群失稳的概率均大于P*，说明继续第二步开采一定会导致矿柱群失稳现象的发生，故在首次回采完成后必须要采取充填和支护等技术措施，消除地压灾害，保障下一步矿柱群安全高效回采。

（1）针对矿柱群失稳自组织现象，基于RG理论建立二维RG模型，采用Weibull分布表征矿柱群失稳规律，推导出矿柱群失稳临界概率表达式，进而计算出矿柱群失稳临界概率和极限临界概率。当矿柱群失稳概率小于临界概率时，矿柱群可实现自稳；当矿柱群失稳概率处于临界概率与极限临界概率之间时，矿柱群开始向失稳趋势演化；当矿柱群失稳概率大于极限临界概率时，矿柱群失稳。

（2）结合长坡—铜坑矿某矿柱群矿柱形态和岩石力学参数，基于多项逻辑回归模型计算出单个矿柱失稳概率，将其代入BN模型，计算出矿柱群失稳概率，通过与临界概率进行对比得出该矿柱群目前处于稳定状态。

（3）通过定义单个矿柱失稳概率对BN模型进行概率反演，并分析单个矿柱对矿柱群的敏感性，选取敏感性最低的待采区域作为首要开采对象。结果表明，当首次回采选择3号待采区时，矿柱群仍处于稳定状态；二次回采时矿柱群失稳概率均大于矿柱群临界失稳概率，故首次回采后不应进行二次开采，同时需采取必要的充填和支护措施。