基于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），在保證結果準確性的前提下提高分析效率。