基于PCA法的岩体震源破坏面产状分析及工程应用

刘胤池 李庶林 周梦婧 唐 超

（厦门大学建筑与土木工程学院，福建厦门361005）

近十多年以来，微震监测技术被广泛用于我国矿山生产、水电工程和隧道工程等安全监测领域，在岩体灾害预警方面取得了一定的应用效果［1-3］。研究工程岩体的震源机制，可以深刻认识工程岩体的细观破裂机理和宏观破坏模式，对于灾害的监测预警具有重要意义。震源机制涉及到岩体破裂源的物理力学机制，是微震监测领域的一个理论核心问题，也是目前微震研究的一个热点。破坏面产状分析是震源机制解中很重要的一部分，常利用矩张量反演或应力场反演来求得破坏面产状［4-5］。

针对求解构造地震破坏面产状的问题，MICHELINI等［6］首次提出了一种快速获取破坏面产状的方法，将主成分分析法运用至地震学中，提出了“局部破裂椭球体”的概念，并根据震源发生的时间顺序及空间分布推断复杂地区微震活动的平均破坏面几何形状以及活动趋势；POSADAS等［7］进一步完善了该方法，补充了空间位置的主成分分析，用这两种方法可以求得地震导致的主破坏面产状以及破坏面的时空演化趋势。上述方法在地震学中得到了较为广泛的应用，得到国内外学者的进一步深入分析、发展和运用［8-11］。

岩石工程中的岩体破裂与地震类似，在宏观破坏面最终形成之前会有大量的微震事件产生，且这些微震事件具有空间和时间上的聚类特征，微震事件在时间以及空间上的分布可以提供关于破坏面空间产状的相关信息。URBANCIC等［12］首次在微震监测中利用主成分分析法快速求解破坏面产状，其结果与原位测量结果相差较小，证明了该方法在矿山定位中的实用性；TRIFU等［13］利用多种方法对矿山的应力状态进行分析，将主成分分析法得到的破坏面构造与P、T轴方向的应力反演结果对比，结果相近；COULSON［14］等运用时间主成分技术分析煤柱屈服区的范围和趋势，通过煤柱中形成断裂结构的事件簇来表征破坏过程；GAMMANS［15］引入平面度与线性度两个概念来描述局部破裂椭球体，并测试了PCA法求解破坏面产状的稳定性。

在地震学中利用PCA法快速求解地震断层面产状信息及演化趋势的应用较多，然而在岩石工程中，仅利用PCA法获取破坏面产状信息，在获取破坏面产状演化趋势方面鲜有研究。为此，本研究将PCA方法应用至柿竹园矿一次典型的特大山体剪切滑移案例中，根据接收到的微震定位数据来获取滑移破坏面的相关信息及演化趋势，以进一步探索PCA方法在岩体剪滑破坏模式中的应用效果。

1 主成分分析原理

主成分分析（PCA）是一种量化数据分布的统计方法，利用正交变换将一系列可能线性相关的变量转换为一组线性不相关的新变量（也称为主成分），利用新变量展示数据的特征［16］。在工程微震定位事件分析中，对定位事件位置应用PCA，可量化微震事件聚簇类的形状（线、平面、球体）和方向的聚集程度。震源定位位置是在一个三维X-Y-Z坐标系中进行标记的，可用一个最佳拟合椭球体来描述定位事件簇的形状。定位事件在空间与时间上的分布会提供有关破坏面的几何形状信息，本研究从空间与时间两个方面来对岩体破坏面进行分析。基于震源空间分布的主成分分析（PCA With a Spatial Window，SPCA）是根据微震信号的空间坐标位置利用PCA获取微震信号的空间分布特征，再进一步获得矿山的主要破坏面的产状；基于震源时间分布的主成分分析（PCA With a Temporal Window，TPCA）将一个覆盖有固定数量事件的时间窗口在按时间顺序排列好的微震事件群上滑动，求得每个窗口代表的椭球体并将方向不同的椭球体进行分离，从而能够获得不同时刻的破坏面产状，利用该方法可以获得破坏面的演化趋势。

1.1 基于震源空间分布的主成分分析

对于空间一定区域内的N个微震定位事件，利用主元分析方法对震源破裂面进行分析，该事件群的协方差矩阵用S表示

该矩阵在几何学上可表示为由震源坐标拟合的空间椭球，其特征向量方向给出了椭球轴的空间方向，特征值为椭球体半轴的长度，该椭球体称为局部破裂椭球体［6，15］。

椭球的椭圆度为e，公式为

当e≥2.5时，椭球体可以看作是微震事件作用下所产生的破裂面；当e＜2.5时，此时的椭球体不能视为平面，而是体积形态，在这种情况下不能根据椭球体的特征来推断破坏面的产状。

对于可以视为平面的椭球体，第1特征向量n1和第2特征向量n2构成的平面即为潜在破坏面，第3特征向量n3为破坏面法向。根据破坏面法向中的3个元素n3（1）、n3（2）和n3（3），可以求得破坏面得出走向φ和倾角δ，公式为［17］

1.2 基于震源时间分布的主成分分析

为了推断岩体潜在破坏面在时间序列上的演化进程，本研究采用基于震源位置时间分布的主成分分析法对监测到的微震定位数据进行分析。

对于N个震源定位事件，将其按时间顺序排列，选定一个长为K的滑动窗口，K为一次性覆盖微震数据的个数，此时每K个事件可以看作一个椭球体，在N个事件上共有p=N-K+1个椭球体。为了减少窗口长度对结果的影响，需选择最佳窗口长度，当Q接近于常数时，此时的K值为最佳窗口长度［6-7］：

将窗口在震源序列上滑动，每个窗口为一个椭球体，椭球体之间主轴的相差角度可进行如下计算：

式中，α1h、β1h分别为第1个椭球体与第h个椭球体主轴相差的角度，（°）；、分别为第h个椭球体的最大特征向量与最小特征向量。

设定分离角度分别为δα与δβ，当满足下列条件时，说明两椭球体主轴方向几乎相同，为分离出主轴不同的椭球体，窗口继续滑动直至式（7）中一个不等式不成立为止，即找出不同主轴方向的椭球体。

当上式不满足时，取该时间段所有椭球体主轴方向向量的平均值作为该时刻椭球体主轴方向，分别用v1(K ,r)和v3(K ,s)表示：

式中，r和s为满足下列不等式的连续椭球体个数。在工程中的微震中，由于事件与事件间相距较近，因此δα与δβ需选择比地震学中更为精确的角度来区分不同的椭球体，本研究两者取值分别为5°和10°。

本研究采用主成分分析法获取破坏面产状及破坏面演化进程信息，计算流程如图1所示。

在两种分析方法中，均按时间顺序将震源坐标进行排列，根据SPCA可以获得局部破裂椭球体的主轴信息，据此可进一步求得破坏面产状；在TPCA中，利用式（5）确定合适的滑动窗口长度K，是为了降低因滑动窗口长度不同对破坏面产状的影响［6］，之后求解每个窗口的局部椭球体信息与SPCA类似，根据式（6）和式（7）可以分离出主轴方向不同的椭球体，从而获得不同时刻的破坏面产状。

2 微震监测案例分析

2.1 工程概况

为了监测采空区、矿柱和上覆厚大岩层的稳定性和安全性，柿竹园有色金属有限责任公司所属的多金属采矿场于2010年引入加拿大ESG的36通道的全数字微震监测系统。该系统全部采用单轴加速度传感器，分别布置在514 m平面、558 m平面和630 m平面内。2012年6月21日该矿进行的一次大爆破落矿过程中，直接诱发了采场上覆岩层（山体）的特大规模滑移型垮塌现象。微震监测系统完整地监测到了这次破坏的孕育、垮塌的全过程［18］。本次大爆破诱发的微震定位事件总共为85个，图2中给出了山体垮塌区域的微震定位事件分布，主要集中于上部悬顶破坏区域的640～755 m高程范围内。这是一次非常典型的、由大爆破诱发的特大山体破坏案例。本研究根据记录到的微震定位数据，利用SPCA法和TPCA法，对矿山破坏面的相关信息进行分析。

2.2 案例分析

根据1.1节和1.2节描述，利用SPCA法与TPCA法能够快速计算出微震活动产生的破坏面产状，并获得破坏面的发展过程。将SPCA法与TPCA法应用于2012年6月21日一次因大爆破产生的特大山体剪切滑移事件中，推断可能产生的潜在破坏面相关信息。本研究选取上部悬顶岩体滑移定位事件进行分析，去除定位误差较大的事件后，最终选取52个定位事件（表1）进行SPCA与TPCA分析。

2.2.1 SPCA法求解破坏面产状

将微震监测系统监测到52个岩体破裂的震源定位事件的坐标按时间顺序进行记录，利用Matlab中的z-score函数对原始空间坐标数据进行标准化处理，结果见表2。

通过对协方差矩阵进行特征值-特征向量分解，并对求出的特征值进行排序。由于前两个主元向量是构成所求潜在破裂平面的基向量，最后一个向量为该平面的法向量，因而使用该法向矢量和位于平面坐标的均值拟合到平面上，计算出可能的破坏面产状。

案例中微震定位数据的协方差矩阵为

求得的特征值从大到小排列分别为λ1=2.484 8,λ2=0.435 6,λ3=0.079 6，其对应的特征向量分别为n1=（-0.610 8，-0.589 8，0.528 3），n2=（0.255 2，0.485 0，0.836 4），n3=（0.749 6，-0.645 7，0.145 7）。

根据1.1节描述，因为λ1/λ3≥ 2.5，即可将该椭球体从几何特征上视为一个平面，该平面即为所求的可能的震源破裂面由n1与n2两轴确定，n3为破坏面法向。

根据式（4）及倾向与走向之间的关系，可求得破坏面的倾向及倾角分别为319.26°与81.62°。在Wulff投影下，用SPCA法求出的主破坏面产状由图3所示；实际测量的破裂面产状倾向为327°和倾角为78°，如图4和图5所示。可见，用SPCA法可以较为快速且精准地获得破坏面产状。

2.2.2 TPCA法求解破坏面产状演化趋势

首先通过Q值的变化对K取值，如图6所示。当K=36时，Q值之后的变化较为缓慢，故选择此时的K值作为滑动窗口长度。

在求得每个窗口所代表的椭球体的主轴方向向量后，利用式（6）和式（7）将不同的椭球体分离出，根据式（8）求得不同椭球体的主轴方向，据此可以得到不同时刻的不同椭球体的破坏面产状，从而获得破坏面的演化趋势。

利用TPCA法获得的破坏面变化信息见表3。

由表3可知：此次滑移事故中，主破坏面倾向从NW316.14°向NW317.47°方向发展，破坏面倾角由87.96°向85.05°发展。

3 结论

本研究分析了基于主成分分析快速获取主破坏面相关信息的方法，该方法分为SPCA和TPCA两部分。根据震源坐标的时间与空间分布，利用SPCA法可以获得局部破裂椭球体的主轴信息，据此可进一步求得破坏面产状；利用TPCA法可以分离出主轴方向不同的椭球体，从而可获得不同时刻破坏面的产状，来判断破坏面的演化趋势。将该方法应用于柿竹园矿山的一次山体滑移事件分析中，快速获取了破坏面的产状及演化趋势，得到了如下结论：

（1）由SPCA法确定了此次山体滑移产生的破坏面产状，其倾向为NW319.26°，倾角为81.62°，与实际测量结果（倾向为NW327°，倾角为78°）较接近，证明了该方法的合理性，此外，SPCA法利用了微震信号的定位信息，因而计算过程较简便。

（2）利用TPCA法分析了破坏面的演化过程，破坏面产状与SPCA法计算结果相近，其倾向变化方向由NW316.14°向NW317.47°方向发展，倾角由 87.96°向85.05°发展，可推断该破坏倾向向北发展，且倾角逐渐减小，可见，TPCA法够推断出破坏面的演化趋势，为灾害预防提供理论支撑。

（3）当定位数据充足时，在求得最佳的窗口长度K后可以进一步考虑最佳步距问题（本研究默认为1），在保证结果准确性的前提下提高分析效率。