基于PCA法的巖體震源破壞面產(chǎn)狀分析及工程應(yīng)用

劉胤池 李庶林 周夢婧 唐 超

（廈門大學(xué)建筑與土木工程學(xué)院，福建廈門361005）

近十多年以來，微震監(jiān)測技術(shù)被廣泛用于我國礦山生產(chǎn)、水電工程和隧道工程等安全監(jiān)測領(lǐng)域，在巖體災(zāi)害預(yù)警方面取得了一定的應(yīng)用效果［1-3］。研究工程巖體的震源機(jī)制，可以深刻認(rèn)識工程巖體的細(xì)觀破裂機(jī)理和宏觀破壞模式，對于災(zāi)害的監(jiān)測預(yù)警具有重要意義。震源機(jī)制涉及到巖體破裂源的物理力學(xué)機(jī)制，是微震監(jiān)測領(lǐng)域的一個理論核心問題，也是目前微震研究的一個熱點(diǎn)。破壞面產(chǎn)狀分析是震源機(jī)制解中很重要的一部分，常利用矩張量反演或應(yīng)力場反演來求得破壞面產(chǎn)狀［4-5］。

針對求解構(gòu)造地震破壞面產(chǎn)狀的問題，MICHELINI等［6］首次提出了一種快速獲取破壞面產(chǎn)狀的方法，將主成分分析法運(yùn)用至地震學(xué)中，提出了“局部破裂橢球體”的概念，并根據(jù)震源發(fā)生的時(shí)間順序及空間分布推斷復(fù)雜地區(qū)微震活動的平均破壞面幾何形狀以及活動趨勢；POSADAS等［7］進(jìn)一步完善了該方法，補(bǔ)充了空間位置的主成分分析，用這兩種方法可以求得地震導(dǎo)致的主破壞面產(chǎn)狀以及破壞面的時(shí)空演化趨勢。上述方法在地震學(xué)中得到了較為廣泛的應(yīng)用，得到國內(nèi)外學(xué)者的進(jìn)一步深入分析、發(fā)展和運(yùn)用［8-11］。

巖石工程中的巖體破裂與地震類似，在宏觀破壞面最終形成之前會有大量的微震事件產(chǎn)生，且這些微震事件具有空間和時(shí)間上的聚類特征，微震事件在時(shí)間以及空間上的分布可以提供關(guān)于破壞面空間產(chǎn)狀的相關(guān)信息。URBANCIC等［12］首次在微震監(jiān)測中利用主成分分析法快速求解破壞面產(chǎn)狀，其結(jié)果與原位測量結(jié)果相差較小，證明了該方法在礦山定位中的實(shí)用性；TRIFU等［13］利用多種方法對礦山的應(yīng)力狀態(tài)進(jìn)行分析，將主成分分析法得到的破壞面構(gòu)造與P、T軸方向的應(yīng)力反演結(jié)果對比，結(jié)果相近；COULSON［14］等運(yùn)用時(shí)間主成分技術(shù)分析煤柱屈服區(qū)的范圍和趨勢，通過煤柱中形成斷裂結(jié)構(gòu)的事件簇來表征破壞過程；GAMMANS［15］引入平面度與線性度兩個概念來描述局部破裂橢球體，并測試了PCA法求解破壞面產(chǎn)狀的穩(wěn)定性。

在地震學(xué)中利用PCA法快速求解地震斷層面產(chǎn)狀信息及演化趨勢的應(yīng)用較多，然而在巖石工程中，僅利用PCA法獲取破壞面產(chǎn)狀信息，在獲取破壞面產(chǎn)狀演化趨勢方面鮮有研究。為此，本研究將PCA方法應(yīng)用至柿竹園礦一次典型的特大山體剪切滑移案例中，根據(jù)接收到的微震定位數(shù)據(jù)來獲取滑移破壞面的相關(guān)信息及演化趨勢，以進(jìn)一步探索PCA方法在巖體剪滑破壞模式中的應(yīng)用效果。

1 主成分分析原理

主成分分析（PCA）是一種量化數(shù)據(jù)分布的統(tǒng)計(jì)方法，利用正交變換將一系列可能線性相關(guān)的變量轉(zhuǎn)換為一組線性不相關(guān)的新變量（也稱為主成分），利用新變量展示數(shù)據(jù)的特征［16］。在工程微震定位事件分析中，對定位事件位置應(yīng)用PCA，可量化微震事件聚簇類的形狀（線、平面、球體）和方向的聚集程度。震源定位位置是在一個三維X-Y-Z坐標(biāo)系中進(jìn)行標(biāo)記的，可用一個最佳擬合橢球體來描述定位事件簇的形狀。定位事件在空間與時(shí)間上的分布會提供有關(guān)破壞面的幾何形狀信息，本研究從空間與時(shí)間兩個方面來對巖體破壞面進(jìn)行分析?；谡鹪纯臻g分布的主成分分析（PCA With a Spatial Window，SPCA）是根據(jù)微震信號的空間坐標(biāo)位置利用PCA獲取微震信號的空間分布特征，再進(jìn)一步獲得礦山的主要破壞面的產(chǎn)狀；基于震源時(shí)間分布的主成分分析（PCA With a Temporal Window，TPCA）將一個覆蓋有固定數(shù)量事件的時(shí)間窗口在按時(shí)間順序排列好的微震事件群上滑動，求得每個窗口代表的橢球體并將方向不同的橢球體進(jìn)行分離，從而能夠獲得不同時(shí)刻的破壞面產(chǎn)狀，利用該方法可以獲得破壞面的演化趨勢。

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

對于空間一定區(qū)域內(nèi)的N個微震定位事件，利用主元分析方法對震源破裂面進(jìn)行分析，該事件群的協(xié)方差矩陣用S表示

該矩陣在幾何學(xué)上可表示為由震源坐標(biāo)擬合的空間橢球，其特征向量方向給出了橢球軸的空間方向，特征值為橢球體半軸的長度，該橢球體稱為局部破裂橢球體［6，15］。

橢球的橢圓度為e，公式為

當(dāng)e≥2.5時(shí)，橢球體可以看作是微震事件作用下所產(chǎn)生的破裂面；當(dāng)e＜2.5時(shí)，此時(shí)的橢球體不能視為平面，而是體積形態(tài)，在這種情況下不能根據(jù)橢球體的特征來推斷破壞面的產(chǎn)狀。

對于可以視為平面的橢球體，第1特征向量n1和第2特征向量n2構(gòu)成的平面即為潛在破壞面，第3特征向量n3為破壞面法向。根據(jù)破壞面法向中的3個元素n3（1）、n3（2）和n3（3），可以求得破壞面得出走向φ和傾角δ，公式為［17］

1.2 基于震源時(shí)間分布的主成分分析

為了推斷巖體潛在破壞面在時(shí)間序列上的演化進(jìn)程，本研究采用基于震源位置時(shí)間分布的主成分分析法對監(jiān)測到的微震定位數(shù)據(jù)進(jìn)行分析。

對于N個震源定位事件，將其按時(shí)間順序排列，選定一個長為K的滑動窗口，K為一次性覆蓋微震數(shù)據(jù)的個數(shù)，此時(shí)每K個事件可以看作一個橢球體，在N個事件上共有p=N-K+1個橢球體。為了減少窗口長度對結(jié)果的影響，需選擇最佳窗口長度，當(dāng)Q接近于常數(shù)時(shí)，此時(shí)的K值為最佳窗口長度［6-7］：

將窗口在震源序列上滑動，每個窗口為一個橢球體，橢球體之間主軸的相差角度可進(jìn)行如下計(jì)算：

式中，α1h、β1h分別為第1個橢球體與第h個橢球體主軸相差的角度，（°）；、分別為第h個橢球體的最大特征向量與最小特征向量。

設(shè)定分離角度分別為δα與δβ，當(dāng)滿足下列條件時(shí)，說明兩橢球體主軸方向幾乎相同，為分離出主軸不同的橢球體，窗口繼續(xù)滑動直至式（7）中一個不等式不成立為止，即找出不同主軸方向的橢球體。

當(dāng)上式不滿足時(shí)，取該時(shí)間段所有橢球體主軸方向向量的平均值作為該時(shí)刻橢球體主軸方向，分別用v1(K ,r)和v3(K ,s)表示：

式中，r和s為滿足下列不等式的連續(xù)橢球體個數(shù)。在工程中的微震中，由于事件與事件間相距較近，因此δα與δβ需選擇比地震學(xué)中更為精確的角度來區(qū)分不同的橢球體，本研究兩者取值分別為5°和10°。

本研究采用主成分分析法獲取破壞面產(chǎn)狀及破壞面演化進(jìn)程信息，計(jì)算流程如圖1所示。

在兩種分析方法中，均按時(shí)間順序?qū)⒄鹪醋鴺?biāo)進(jìn)行排列，根據(jù)SPCA可以獲得局部破裂橢球體的主軸信息，據(jù)此可進(jìn)一步求得破壞面產(chǎn)狀；在TPCA中，利用式（5）確定合適的滑動窗口長度K，是為了降低因滑動窗口長度不同對破壞面產(chǎn)狀的影響［6］，之后求解每個窗口的局部橢球體信息與SPCA類似，根據(jù)式（6）和式（7）可以分離出主軸方向不同的橢球體，從而獲得不同時(shí)刻的破壞面產(chǎn)狀。

2 微震監(jiān)測案例分析

2.1 工程概況

為了監(jiān)測采空區(qū)、礦柱和上覆厚大巖層的穩(wěn)定性和安全性，柿竹園有色金屬有限責(zé)任公司所屬的多金屬采礦場于2010年引入加拿大ESG的36通道的全數(shù)字微震監(jiān)測系統(tǒng)。該系統(tǒng)全部采用單軸加速度傳感器，分別布置在514 m平面、558 m平面和630 m平面內(nèi)。2012年6月21日該礦進(jìn)行的一次大爆破落礦過程中，直接誘發(fā)了采場上覆巖層（山體）的特大規(guī)?；菩涂逅F(xiàn)象。微震監(jiān)測系統(tǒng)完整地監(jiān)測到了這次破壞的孕育、垮塌的全過程［18］。本次大爆破誘發(fā)的微震定位事件總共為85個，圖2中給出了山體垮塌區(qū)域的微震定位事件分布，主要集中于上部懸頂破壞區(qū)域的640～755 m高程范圍內(nèi)。這是一次非常典型的、由大爆破誘發(fā)的特大山體破壞案例。本研究根據(jù)記錄到的微震定位數(shù)據(jù)，利用SPCA法和TPCA法，對礦山破壞面的相關(guān)信息進(jìn)行分析。

2.2 案例分析

根據(jù)1.1節(jié)和1.2節(jié)描述，利用SPCA法與TPCA法能夠快速計(jì)算出微震活動產(chǎn)生的破壞面產(chǎn)狀，并獲得破壞面的發(fā)展過程。將SPCA法與TPCA法應(yīng)用于2012年6月21日一次因大爆破產(chǎn)生的特大山體剪切滑移事件中，推斷可能產(chǎn)生的潛在破壞面相關(guān)信息。本研究選取上部懸頂巖體滑移定位事件進(jìn)行分析，去除定位誤差較大的事件后，最終選取52個定位事件（表1）進(jìn)行SPCA與TPCA分析。

2.2.1 SPCA法求解破壞面產(chǎn)狀

將微震監(jiān)測系統(tǒng)監(jiān)測到52個巖體破裂的震源定位事件的坐標(biāo)按時(shí)間順序進(jìn)行記錄，利用Matlab中的z-score函數(shù)對原始空間坐標(biāo)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，結(jié)果見表2。

通過對協(xié)方差矩陣進(jìn)行特征值-特征向量分解，并對求出的特征值進(jìn)行排序。由于前兩個主元向量是構(gòu)成所求潛在破裂平面的基向量，最后一個向量為該平面的法向量，因而使用該法向矢量和位于平面坐標(biāo)的均值擬合到平面上，計(jì)算出可能的破壞面產(chǎn)狀。

案例中微震定位數(shù)據(jù)的協(xié)方差矩陣為

求得的特征值從大到小排列分別為λ1=2.484 8,λ2=0.435 6,λ3=0.079 6，其對應(yīng)的特征向量分別為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）。

根據(jù)1.1節(jié)描述，因?yàn)棣?/λ3≥ 2.5，即可將該橢球體從幾何特征上視為一個平面，該平面即為所求的可能的震源破裂面由n1與n2兩軸確定，n3為破壞面法向。

根據(jù)式（4）及傾向與走向之間的關(guān)系，可求得破壞面的傾向及傾角分別為319.26°與81.62°。在Wulff投影下，用SPCA法求出的主破壞面產(chǎn)狀由圖3所示；實(shí)際測量的破裂面產(chǎn)狀傾向?yàn)?27°和傾角為78°，如圖4和圖5所示。可見，用SPCA法可以較為快速且精準(zhǔn)地獲得破壞面產(chǎn)狀。

2.2.2 TPCA法求解破壞面產(chǎn)狀演化趨勢

首先通過Q值的變化對K取值，如圖6所示。當(dāng)K=36時(shí)，Q值之后的變化較為緩慢，故選擇此時(shí)的K值作為滑動窗口長度。

在求得每個窗口所代表的橢球體的主軸方向向量后，利用式（6）和式（7）將不同的橢球體分離出，根據(jù)式（8）求得不同橢球體的主軸方向，據(jù)此可以得到不同時(shí)刻的不同橢球體的破壞面產(chǎn)狀，從而獲得破壞面的演化趨勢。

利用TPCA法獲得的破壞面變化信息見表3。

由表3可知：此次滑移事故中，主破壞面傾向從NW316.14°向NW317.47°方向發(fā)展，破壞面傾角由87.96°向85.05°發(fā)展。

3 結(jié)論

本研究分析了基于主成分分析快速獲取主破壞面相關(guān)信息的方法，該方法分為SPCA和TPCA兩部分。根據(jù)震源坐標(biāo)的時(shí)間與空間分布，利用SPCA法可以獲得局部破裂橢球體的主軸信息，據(jù)此可進(jìn)一步求得破壞面產(chǎn)狀；利用TPCA法可以分離出主軸方向不同的橢球體，從而可獲得不同時(shí)刻破壞面的產(chǎn)狀，來判斷破壞面的演化趨勢。將該方法應(yīng)用于柿竹園礦山的一次山體滑移事件分析中，快速獲取了破壞面的產(chǎn)狀及演化趨勢，得到了如下結(jié)論：

（1）由SPCA法確定了此次山體滑移產(chǎn)生的破壞面產(chǎn)狀，其傾向?yàn)镹W319.26°，傾角為81.62°，與實(shí)際測量結(jié)果（傾向?yàn)镹W327°，傾角為78°）較接近，證明了該方法的合理性，此外，SPCA法利用了微震信號的定位信息，因而計(jì)算過程較簡便。

（2）利用TPCA法分析了破壞面的演化過程，破壞面產(chǎn)狀與SPCA法計(jì)算結(jié)果相近，其傾向變化方向由NW316.14°向NW317.47°方向發(fā)展，傾角由 87.96°向85.05°發(fā)展，可推斷該破壞傾向向北發(fā)展，且傾角逐漸減小，可見，TPCA法夠推斷出破壞面的演化趨勢，為災(zāi)害預(yù)防提供理論支撐。

（3）當(dāng)定位數(shù)據(jù)充足時(shí)，在求得最佳的窗口長度K后可以進(jìn)一步考慮最佳步距問題（本研究默認(rèn)為1），在保證結(jié)果準(zhǔn)確性的前提下提高分析效率。