基于尖点突变模型与D-S证据融合理念的地下矿山岩体失稳预警方法

罗周全，李艳艳，秦亚光，文 磊

(中南大学资源与安全工程学院，湖南 长沙 410083)

0 引言

岩体失稳是指在应力作用下岩体内部微裂隙不断发育形成断裂面，断裂面贯通后岩体瞬间破坏脱落的过程[1-3]。深部开采岩体失稳灾害发生频率高、破坏性强，开展岩体失稳预警研究对于保障矿山经济效益和人员安全意义重大[3-6]。目前对岩体失稳灾害的预警多借助于微震监测、声发射检测、红外热像、微重力检测等手段[7-11]。其中，利用微震监测技术采集并分析岩体变形破坏时释放的弹性波，提取特征参数对岩体失稳进行预警的方法被广泛应用。李瑞等[12]对大量微震监测数据进行分析，发现视体积、能量指数等参数在岩体失稳之前会呈现异常变化趋势，利用这些参数的特性可以进行岩体失稳预警。王俊光等[13]通过对微震统计数据进行时序分析，揭示了微震能量、频次、振动速度和幅值等参数在冲击地压发生前的活动规律，对矿山地质灾害预警具有一定的指导意义。张楚旋等[14]结合沙坝矿实际微震监测参数，利用b值、能量指数、施密特数和累积视体积建立了岩体失稳的智能预报模式。

上述研究明确了一系列对岩体稳定状态敏感的微震特征参数，但是利用其进行岩体失稳预警时多偏向于定性分析，预警结果的准确性和可靠性均有待进一步提高。尖点突变理论能够有效利用岩体失稳灾害具有的突发性特点，定量分析岩体的稳定状态并对失稳灾害进行预警。张钦礼等[15]通过尖点突变理论对建立的简化力学模型进行分析，推导出采场顶板-矿柱系统失稳的临界条件。付成华等[16]应用能量、熵、位移模等参数建立了5种判据对围岩稳定性进行分析，从不同方面对岩体失稳进行预警。刘新荣等[17]结合有限元方法分析边坡强度折减过程中能量的变化规律，建立了适用于边坡失稳预警的能量突变判据。综上，利用尖点突变预警模型能实现对岩体稳定性进行定量分析，且模型具有自适应鲁棒性，缺点是突变判据单一时预警准确性难以保障，突变判据较多时又易发生冲突，以致模型的预警结果准确性难以保障。

D-S证据融合理论是一种可以有效解决证据冲突的多参数评价模型。袁四化等[18]依据时空观和借助D-S证据融合理论建立了稳定性评价多指标体系和多信息融合评价模型，为矿山采空区稳定性评价提供了借鉴。LUO Z Q等[19]从岩体失稳的内部机理和外部表现入手，通过将多场耦合分析得到损伤指标并与微震特征参数指标作证据融合建立起岩体失稳预警模型。罗小燕等[20]采集单轴压缩实验各阶段的声发射特征参数，运用敏感度特征提取法等方法选定证据体，将各证据体融合后获得各证据体于各阶段的信任度权重，建立加权证据体融合的失稳预警模型。综上，D-S证据融合理论可以用于岩体失稳灾害的预警，但是D-S证据融合理论在运用时仍存在鲁棒性较差、模型抗干扰能力弱的缺点。

鉴于尖点突变理论与D-S证据融合理论两种方法在运用时呈现出各有不足却优缺点互补的现象。笔者提出一种将尖点突变分析与D-S证据融合预警综合运用的新的预警方法，力求提高岩体失稳预警结果的准确性。本文依托广东某地下金属矿山为研究背景，基于由尖点突变理论和D-S证据融合理论形成的综合预警方法对-550 m采矿中段至-650 m采矿中段的岩体稳定性进行分析与预警，以期对矿山岩体失稳灾害的有效预警提供新的研究思路。

1 基本原理

1.1 尖点突变模型

尖点突变模型是突变理论中的一种初等突变模型，它可以忽略系统尚未可知的内部作用而直接处理飞跃式的不连续变化。主要用于研究复杂的非线性变化过程，判断系统突变发生的条件，评价系统的稳定性[21-24]。该方法在岩体稳定性评价及失稳灾害的预警中广泛应用。岩体在破裂时会通过微地震活动释放能量[25]，因此微震事件的活动性参数与岩体状态变化之间具有密切相关性，岩石状态的突变也对应着微震活动性参数的突变。能量指数EI作为微震活动性参数之一，亦具有突变性，符合尖点突变分析的使用条件。

根据微震能量指数EI的监测序列获得能量指数的岩体失稳突变函数，并推导出尖点突变模型的突变特征值Δ的表达式为：

Δ=8u3+27v2

(1)

式中，u和v为基于能量指数的岩体失稳尖点突变势函数的系数参数。详细推导过程见孙建等[25]的研究。Δ值的大小，代表岩石处于不同的稳定状态，即：

(1)Δ>0时，岩体只对应一种状态，即稳定状态。代表此时岩体承受的应力未达到峰值，微裂隙正在发育中，岩体变形还处于连续变化状态，稳定性相对较强；

(2)Δ=0时，岩体处于临界状态，稍加干扰就可能从稳定状态突跳到不稳定状态，产生失稳风险，岩体稳定性欠佳。

(3)Δ<0时，岩体发生突变的可能性较大，岩体稳定性较差。

利用尖点突变模型进行岩体失稳预警时，判据Δ<0且Δ值越小失稳风险越大，在一段时间内Δ<0且Δ值最小的时间内岩体失稳发生的可能性最大。

1.2 D-S证据融合理论

D-S证据理论是于1967年被提出的一种不精确推理理论。该理论是在证据理论的基础上利用一定的证据组合规则将各个证据体进行融合，并根据相应规则做出最终决策[26]。作为一种处理不确定信息的方法，D-S证据融合理论近年来被广泛应用于边坡、岩体稳定性评价方面[27]。因此，本文采用D-S融合理论对岩体稳定性评价是完全可行的。

日最大岩层振动速度(Vmax)、日视体积增量(δVA)与施密特指数(Sc)是微震监测常用的参数，也被广泛应用到矿岩活动性评价中[28]。本文选取这三个参数作为证据体进行融合分析。利用微震监测结果获取三个参数的统计数据，依据D-S证据融合理论建立多参量数学预警模型，以科学合理地评价对象区域的岩体稳定性。

(1)风险辨识框架设定

设定证据集Ω={ Λ1，Λ2，Λ3}，其中，Λ1代表Vmax，Λ2代表δVA，Λ3代表Sc。

根据岩体失稳风险等级分为3个级别，各个证据体的不同取值对应不同的岩体失稳风险等级，即：

Ⅰ级风险：表示监测对象区域岩体属于低失稳风险，岩体失稳发生的概率非常小；

Ⅱ级风险：表示监测对象区域岩体属于中等失稳风险，岩体失稳发生的概率较小；

Ⅲ级风险：表示监测对象区域岩体属于高失稳风险，岩体失稳发生的概率较大。

设αn，βn为证据Λn所对应的划分Ⅰ级风险、Ⅱ级风险和Ⅲ级风险的临界值(图1)。

图1 依据各参数的风险等级划分Fig.1 Risk level division according to the parameter

证据体Λn对应于不同风险等级的取值叫做不同证据体对应于各风险等级的基本可信度(BPA)，符号记做Λnk。

Λnk(n=1,2,3;k=Ⅰ,Ⅱ,Ⅲ)

(2)

在式(2)中，n=1时，对应的证据体为日最大岩层振动速度(Vmax)；n=2时，对应的证据体为日视体积增量(δVA)；n=3时，对应的证据体为施密特指数(Sc)；而(k=Ⅰ,Ⅱ,Ⅲ)则分别对应3种不同等级风险。Λnk的求值公式如式(3)～式(5)。

(3)

(4)

(5)

根据Λnk的求值公式(3)～(5)求出的各证据体对应各等级风险的基本可信度取值，具体结果的表达见表1。

表1 基本可信度赋值Table 1 Basic possibility assignment

(2)证据融合

各个证据体的基本可信度赋值获得后，两个证据相互之间的相似系数ρ可以通过计算得到。ρ的取值区间为[0,1]，ρ值越大，表明两个证据相互之间相互验证性越强，对评价结果的支持度越高，可靠度越大。

证据Λ1和Λ2之间的相似系数计算如下：

(6)

将各证据之间的相似系数的计算结果以矩阵形式阵列出来，称之为相似矩阵，即：

(7)

根据相似矩阵计算单个证据体对最终评价结果的支持度U(Λn)，然后进行归一化处理，使得每个证据体所对应所有等级风险值之和为1，即可得到各个证据体的可信度C(Λn)，计算公式如下：

(8)

(9)

为了更准确地评价岩体状态，将可信度作为证据Λn的权重，结合证据体Ω中各个证据体的基本可信度值解得各个证据对结果的加权合成后的风险概率值，即融合风险概率Λ(k)：

(10)

当Λ(k)值最大时，k所对应的风险等级即为风险等级评价结果。利用D-S证据融合模型进行岩体失稳预警时，若出现失稳等级不断上升并出现1次以上高等失稳风险，则可判断高等失稳风险出现后的一段时间岩体失稳发生的可能性较大。

1.3 综合预警

风险矩阵法是一种常用于风险评估领域的风险可视化方法，其核心是通过一个二维表格实现对风险等级的综合评估。上述尖点突变分析模型和D-S证据融合模型两种模型的风险等级评价结果均在一定程度上反映了对象区域的岩体失稳风险性，但是单依靠其中任一结果都无法保证岩体失稳预警的准确性。为了解决这个问题，进一步提高岩体失稳预警的准确性，以尖点突变分析模型和D-S证据融合模型的评价结果作为两种评估依据，运用风险矩阵对两种结果进行融合，获得综合性预警模型的风险等级表达。基于风险矩阵法的岩体失稳风险等级综合评价矩阵见表2。

表2 综合预警模式Table 2 Comprehensive warning mode

1级：无风险，表示岩体处于稳定的状态，是岩体失稳的早期孕育期，微震活动性较弱。微震活动参数处于原有微裂隙被压缩闭合阶段，但是未产生塑性变形。阶段应对措施：正常生产。

2级：低等风险，岩体处于较稳定的状态，是岩体失稳的缓慢发育期，微震活动性增强。阶段应对措施：正常生产，适当加强排查。

3级：中等风险，岩体处于欠稳定状态，是岩体失稳的快速发育期，微震活动性较强。阶段应对措施：正常生产，加强排查，必要时加强支护与防控。

4级：高等风险，岩体处于稳定与不稳定的交叉临界状态，是岩体失稳的发育峰值期，是岩体能量积聚达到承受极限的状态。阶段应对措施：减缓生产，全范围排查，必要时加强支护与防控。

5级：极高风险，岩体处于不稳定状态，是岩体失稳的能量释放期，极易发生失稳灾害。阶段应对措施：对应区域立即停产，同时加强支护与防控。

利用综合预警模型进行岩体失稳预警时，若出现失稳等级不断上升并出现1次以上4级(高等失稳风险)，则判断高等失稳风险出现后的一段时间岩体失稳发生的可能性较大；5级(极高等失稳风险)一旦出现，就要发出岩体失稳灾害即将发生的警告，呼吁矿方立即采取相应措施。

2 应用实例

2.1 工程概况

广东某地下金属矿山开采已由浅部转入深部，其深部采区的主矿体赋存于F3断层的下盘，埋深-580～-880 m 之间，矿体以块状产出，上、下盘围岩均为灰岩。与浅部相比，深部矿岩较脆，应力集中加剧，地压活动增多，岩体失稳风险上升。且F3断层的存在，加剧了矿区井下地压活动特征的复杂性，增加了地压灾害事故产生的概率。为了确保深部矿体安全、高效地开采深部矿体。该矿山构建微震地压监测系统，对矿山地压活动开展实时监测，并对地压灾害事件进行预警预报。

根据矿山实际生产情况，整个矿山总共布置44个通道微震监测系统，配置32个微震检波器，覆盖10个中段。考虑到目前正开采-550～-650 m之间的3个中段，总段高100 m，控制矿体长度800 m，控制宽度约200 m，因此重点监测该区域。监测该区域的微震检波器主要安装在-550 m和-650 m两个中段(图2)。每个中段布置5个单向检波器(U)和一个三向检波器(T)，共12个检波器，能完全实现该中段的监测。

图2 检波器阵列分布图Fig.2 Geophone array distribution

图2中，红色原点代表该中段工作站所在位置；三角形符号所位置即为检波器布置位置，三角符号对应编号代表检波器的位置和属性信息，如“-550S21U”即代表该检波器为安装在-550m中段，编号为S21的单向检波器；黑色字体如“N10穿脉”、“S1穿脉”代表对应穿脉的在中段中的编号。

2.2 数据获取

系统建成后5个月的统计结果显示，地压活动性最强的区域为-550 m中段～-650 m中段，同时该区域于2017年7月26日发生了具有一定规模的冒顶片帮事件，并且微震震级达1.4级。因此，选取该中段为研究对象，并将“7·26”事件作为已知事件进行尖点突变分析与D-S证据融合分析，以反向验证本文所建立模型的有效性。依据该思路，选择6月31日～7月27日时间段内能量指数监测数据做尖点突变分析，统计结果如图3所示。选择6月31日～7月27日时间段内最大岩层振动速度、日视体积增量与施密特指数的监测数据做D-S证据融合分析，统计结果见表3。

图3 各组能量指数对应序列拟合曲线图Fig.3 Fitting curve corresponding to energy index

2.2.1尖点突变公式拟合

以能量指数为状态变量进行尖点突变分析时,首先要选定用于拟合突变函数的原始数据序列。选定规则为：

(1)以微震系统10天的监测值形成原始数据序列，且序列随日期更新。

(2)更新时剔除前一组序列中最早日期的监测值，添加最新日期的监测值至序列末尾形成新一组序列。

获取2017年6月30日～2017年7月27日期间该矿山微震能量指数的监测数据，并按照上述序列生成规则形成初始序列见表3。由表3可知，能量指数在所选时间段内的监测值所处区间为[-0.25，0.3]，且存在前期随时间逐渐增大，后期随时间逐渐减小的规律，即监测数据序列存在尖点，符合尖点突变分析方法的使用条件。

表3 能量指数序列Table 3 Sequences of energy index

根据表3各组数据的序列次序变化绘制折线图，并对其进行四次多项式拟合，得到如图3所示的一系列拟合曲线。从图3可以看出，7月21日之前各组序列均呈持续上升趋势，7月22日开始序列开始先上升后下降，出现尖点，进一步说明能量指数具有进行尖点突变分析的可行性。

从曲线拟合的报表中可以得到拟合曲线对应的函数表达式以及拟合度R2。从R2的值均大于0.9可知，曲线的拟合效果较优，满足进行尖点突变分析的需要。根据上述分析，采用尖点突变方法，拟合函数表达式进行计算，求出突变方程控制变量u、v的值见表4。

将表4中的控制变量u、v按照式(10)计算突变特征值Δ(图4)。由图4可知，7月9日～7月18日期间，突变特征值Δ持续大于0，说明该监测时间段-600 m中段岩体正处于稳定状态。7月18日～7月27日期间，Δ值多次出现负值，且负值最小值出现在7月22日。由此判断，7月22日发生岩体失稳的可能性最大，但7月22日并未发生岩体失稳事件。这说明尖点突变模型对岩体失稳的预警效果较差。

图4 能量指数突变分析结果Fig.4 Analysis results of energy index based on Cusp catastrophe theory

表4 拟合结果Table 4 Fitting results

2.2.2D-S证据融合预警

统计矿山2017年1月～2017年6月的微震监测数据，分析证据体Vmax、δVA、Sc每日监测值与日最大震级之间的对应关系，确定划分失稳风险等级的临界值αn，βn，具体取值见表5。

以表5中的αn，βn值为参考，将各个证据体不同时期的数据代入公式3～5计算Λnk，获得的Λnk的值即为各证据体不同时期对应的基本可信度取值。结合图5中各证据体于6月30日～7月27日期间的监测数据，计算得到各个证据体的基本可信度并作出其随时间变化的曲线(图6)。由图6知，3个证据体对各级风险的基本可信度取值不同，即不同证据体对风险等级的评估结果存在差别。只有将所有孤立的评估结果综合起来进行分析才能实现直观的、准确的、系统的风险预警。图7所示即为将3个证据体的基本可信度依据式(6)～(10)综合计算得到的融合结果。

图5 各证据体演化规律图Fig.5 Variation of the evidences

图6 证据体基本可信度变化规律Fig.6 Variation of basic possibility of evidences

表5 不同证据体划分各级失稳风险等级的临界值Table 5 Critical values based on different evidence for risk levels of rock mass failure

图7所示为将证据体Vmax、δVA、Sc的等级评价结果融合后得到的各级别风险的融合概率。落入灰色方框区域内的点对应的风险等级即判定为当天的失稳风险等级。从图7中可以得到如下判断：6月30日～7月20日岩体一直处于低失稳风险状态；7月21日失稳风险上升至中等后，于7月22日转换为低风险，又于7月23日～7月25日之间上升至中等风险，据此无法判断7月25日之后是否会发生岩体失稳灾害。进一步验证了D-S证据融合方法鲁棒性较差，无法保证岩体失稳灾害的预警效果。

图7 D-S证据融合预警结果Fig.7 Warning results based on D-S fusion evidence

2.2.3综合预警

基于尖点突变分析与D-S证据融合的风险评价结果，按照表2中的风险矩阵法对岩体失稳风险等级进行综合预警，预警结果与当日最大微震等级见表6。

表6 综合预警结果与实况Table 6 Comprehensive warning results and theactual situations

从表6中可以看出：综合预警模型的预警结果显示7月9日～7月14日以及7月16日～7月17日失稳风险等级为1(无风险)，仅7月15日失稳风险等级为2(低风险)，表明该段时间内监测中段的岩体稳定性较好，无失稳风险；7月18日～7月23日时间段内，失稳风险等级在2级(低风险)和3级(中等风险)之间转换3次，表明该段时间内岩体稳定性变化频繁，且失稳风险等级整体上升；7月24日、7月25日两天时间内岩体失稳风险均上升至4级(高等风险)，已经达到岩体失稳的发育峰值期，表明岩体的稳定状态可能已经比较脆弱，因此判断7月25日之后监测区域发生岩体失稳灾害的可能性极大。

实际情况为：2017年7月26日下午16∶50分左右，矿山监测到一个里氏震级1.4级的特大震级微震事件。矿方在图2中N5、N6穿脉附近发现了具有一定规模的冒顶和片帮现象，巷道顶板有新鲜断面和剥落痕迹(图8)。此外，表5中7月27日失稳风险等级降低至高等风险，这与大震级事件发生后，对应区域的应力和能量得到了部分释放，岩体稳定性上升的现象相一致。

图8 破坏现场照片Fig.8 Photos of rock mass failure

3 结论

井下矿山岩体失稳的准确预警对于矿山安全生产具有重要的意义。本文针对尖点突变分析、D-S证据融合理论两种方法综合预警岩体失稳灾害，得到如下结论：

(1)分别基于尖点突变分析和D-S证据融合理论对岩体失稳灾害进行预警分析，结果表明两种方法的预警都有一定的缺陷。

(2)针对尖点突变分析预警指标过于单一和D-S证据融合理论普适性较弱的缺陷，在利用实时微震监测数据拟合获得突变分析结果的基础上，结合D-S证据融合方法构建了关于岩体失稳风险的综合预警模型。

(3)综合预警模型实现了岩体失稳风险等级的定量化表达，重新设计了失稳风险等级。

(4)结合微震监测数据对综合预警模型的预警效果进行验证，结果表明：综合预警模型有效克服了尖点突变分析和D-S证据融合理论模型单独应用时存在的缺陷，预警效果准确度更高且时效性也更好。

综上，尖点突变与D-S证据融合的失稳综合预警模型为岩体稳定性评价与失稳灾害的预测预警提供了新的技术支撑。