陶 慧，李 瑩，馬小平

（1.河南理工大學電氣工程與自動化學院，河南 焦作 454000；

2.中國礦業大學信息與電氣工程學院，江蘇徐州 2211162；3.焦作大學計算機學院，河南 焦作 454000）

基于多變量監測時序的沖擊地壓復雜性分析

陶 慧1，2，李 瑩3，馬小平2

（1.河南理工大學電氣工程與自動化學院，河南 焦作 454000；

2.中國礦業大學信息與電氣工程學院，江蘇徐州 2211162；3.焦作大學計算機學院，河南 焦作 454000）

論文基于多變量時間序列相空間重構來計算數據的關聯維數，以研究沖擊地壓監測數據的復雜程度。考慮到沖擊地壓監測數據含有噪聲而且長度有限，對傳統G-P算法進行了擴展改進，給出了改進算法求解多變量時間序列的關聯維數的原理，并用于Lorenz混沌系統檢驗了改進算法的有效性。收集了不同沖擊情況下多種監測類型的沖擊地壓時間序列數據，用改進G-P算法求解這些監測數據的關聯維數值。研究結果表明：沖擊地壓監測數據具有混沌特性，而且數據關聯維數越大，復雜程度越高，對應礦井的沖擊破壞性越強。這為基于混沌理論預測沖擊危險性提供了新方法和依據。

沖擊地壓；多變量時間序列；關聯維數；混沌特性

沖擊地壓是礦井生產中最嚴重的自然災害之一，其破壞可比擬為一類非線性微分方程的倍周期分叉而出現的混沌運動現象，謝和平、尹光志等已證實其動力學行為具有混沌特性［1－2］。沖擊破壞過程十分復雜，很難建立精確的數學模型，監測方法有微震、電磁輻射、鉆屑法等［3－5］，可以從現場獲取多種沖擊地壓監測時間序列，通過觀測時間序列來分析其復雜混沌特性具有重要意義。描述系統混沌特性的幾何不變量主要指關聯維數和Lyapunov指數。其中關聯維數d2是系統復雜性的度量，大于關聯維數的下一個整數是刻畫系統所需的獨立變量的個數，為從時間序列恢復原始系統確定了一個框架。大量學者對沖擊地壓監測時間序列的關聯維數d2進行了分析計算，以描述數據的復雜度。如文獻［6］對不同尺度的電磁輻射數據的混沌特性進行了分析，文獻［7］研究了微震時間序列的關聯維數等混沌特性。

但目前對沖擊地壓監測數據的混沌特性分析都基于單變量時間序列。根據F.Takens的嵌入理論，在數據無限長、無噪聲的理想條件下，單變量時間序列足以重構原動力系統［8］，但實際數據一般并不能滿足理想條件，實踐中并不能保證任何單變量時間序列都能很好的實現重構。另外，由于沖擊地壓演化過程復雜，用單一指標進行識別本身就存在很多不足［9］，多變量時序重構包含更加豐富的系統信息，可以彌補數據長度的不足，克服噪聲的影響，因此基于多變量時序重構的關聯維數計算結果更加真實準確［10］。

最常用計算關聯維數d2的方法是G-P算法［11］，一般只針對單變量時間序列相空間重構。本文將對G-P算法進行擴展改進后求解多變量時序重構的關聯維數d2；接下來將改進的G-P算法用于標準混沌系統以說明其有效性，最后將其用于求解沖擊地壓監測時序的關聯維數d2，以描述數據的復雜程度。

1 改進G-P算法計算多變量時序的關聯維數

關聯維數是系統復雜性的一個很好的度量，混沌系統的關聯維數為分數維數，一般認為，大于關聯維數的下一個整數是刻畫系統所需的獨立變量的個數，這為從時間序列恢復原始系統確定了一個框架［8，10］。對單變量時間序列，最常用的關聯維數計算方法是 G-P算法［11］，本文將 G-P算法進行擴展以求解多變量時序相空間重構的關聯維數d2。

對長度為N的有限時間序列，重構的狀態變量個數Nm＝N－N0＋1，是有限的，而關聯積分的值只有在N充分大的情況下接近理論值。常用的計算關聯維數方法是讓嵌入維數m變化，來觀測雙對數lnCN（r）-lnr圖的斜率。當m＜d2時，由于重構相空間的維數不夠高，一般不能分辨吸引子結構，此時lnCN（r）-lnr圖的斜率近似為嵌入維數；當m增加到m＞d2時，吸引子的分辨能力改善，斜率達到一個平穩時期，在雙對數lnCN（r）-lnr圖中找到平行直線區域，這些直線的斜率與m的依賴性消失，該斜率即為關聯維數d［12］。

2

考慮到多變量重構時，嵌入維數m增加對應于多種mi的變化，本文將依次增加各個重構變量的嵌入維數，讓嵌入維數m在一定的范圍內變化，來確定雙對數lnCN（r）-lnr圖的線性部分，再通過斜率計算對應的關聯維數。在實際應用中，確定了線性部分后，為了消除計算誤差，一般采用最小二乘法來擬合雙對數lnCN（r）-lnr圖的線性部分的斜率。

取雙對數線的線性區段內的K個點r1，r2，…rk，由最小二乘法得到式（6）。

2 改進G-P算法應用實例

下面先通過一個具體的例子來說明改進的G－P算法求解多變量時間序列關聯維數的計算效果。對文獻［12］的Lorenz系統，取數據長度為3000并進行歸一化，分別利用第一個變量或者前兩個變量進行相空間重構，并在原始數據的基礎上添加5%的噪聲，得到的雙對數lnCN（r）-lnr圖見圖1。

在圖1中，單變量重構時，重構參數m＝2∶1∶10，τ＝3。圖1（a）表示采用無噪聲x1進行重構的結果；圖1（b）為附加5%噪聲的x1進行重構的結果；雙變量重構時，令τ1＝τ2＝3，從m1＝2，m2＝1依次增加m1和m2使m＝3∶2∶19；圖1（c）表示采用無噪聲的x1和x2進行雙變量重構的結果；圖1（d）表示采用附加5%噪聲的x1和無噪聲的x2進行雙變量重構的結果。

圖1 求解Lorenz的關聯維數的lnCN（r）-lnr圖

分析圖1可以看出，即使數據含有一定的噪聲，在相同的嵌入維數下，雙對數曲線的變化也不大，能比較容易找到雙對數曲線的線性部分。在m＝5時取20個點進行最小二乘擬合，求得四種重構方式下對應的關聯維數如表1所示。根據參考文獻［13］，將d2＝2.04作為Lorenz系統的理論關聯維數值。由于本文數據長度較短，采用單變量重構求得的嵌入維數有一定誤差，但采用多變量時間序列重構時計算結果非常精確，但本文所用數據的長度比參考文獻［13］中要短得多；另外，噪聲的存在使求得的關聯維數d2的值略有增加，但增加不大。

表1 Lorenz系統的關聯維數

總之，采用改進的G-P算法估計多變量時序的關聯維數，對噪聲干擾的魯棒性好，對數據長度要求比單變量時序低，能夠方便且準確地求出基于多變量時序重構的混沌系統的關聯維數，描述系統的復雜程度。

3 沖擊地壓數據來源

為了全面分析沖擊地壓監測數據的混沌特性，本文收集了不同沖擊情況下的多種監測數據。

3.1 微震時間序列數據

煤巖體發生破裂、滑動過程中，會產生一定頻率的向周圍傳播的地震波，利用微震監測系統，可確定各種震動參數。微震事件與沖擊地壓的孕育和發生密不可分［3］。微震能量釋放值能夠表現出微震活動的階段性和微震系統整體特征指標，因此將一段時間區間的累計微震能量指標作為一個微震監測變量。另外，微震時間分布具有層次結構和自相似性，因此將一段時間區間的微震頻次指標，作為第二個微震監測變量。根據微震監控系統我們可以獲得累計能量、累計頻次2個微震監測時間序列。需要說明的是，由于微震能量的范圍從101到108，變化幅度較大，為了更好的反映其變化規律，對其做了取對數處理。

本文微震監控時間序列來源于3個礦井，分別為平煤十一礦，義馬千秋煤礦、義馬躍進煤礦。平煤十一礦由于存在多個孤島工作面，局部應力高度集中，從2008年至今發生40余次沖擊地壓事故，但微震能量都不高，沖擊破壞性很小，本文收集了該礦2009／10／01至2012／12／30期間每天微震累計能量和累計頻次兩個時間序列作為第一組微震數據。義馬千秋煤礦主采煤層為侏羅系二煤，具有中等沖擊傾向性，基本頂為厚度205m的巨厚礫巖層，該礦井的沖擊地壓威脅十分嚴重，共計發生破壞性沖擊地壓20多次，其中破壞巷道超過60m的沖擊次數占60%以上，相比于平煤十一礦，沖擊破壞性更強。本文整理得到的2010／04／06至2012／09／23每天累計能量和頻次微震時間序列作為第二組微震數據。義馬躍進煤礦與千秋煤礦情況相似，至2004年首次發生沖擊地壓以來，該礦已發生沖擊地壓30多次，整理的2010／06／01到2012／9／30日的每天累計能量、累計頻次兩個時間序列作為第3組微震數據。

3.2 電磁輻射時間序列

煤巖體受載發生破裂，因裂縫的形成及擴展過程中電荷的遷移及煤巖顆粒的摩擦而產生電磁輻射。電磁輻射的變化反映了煤巖破壞發生發展的過程，電磁輻射強度（幅值）反映了煤巖體受載及變形破裂強度，脈沖數反映了煤體應力集中程度［4］。一般采用KBD5型便攜式電磁輻射監測儀，在工作面及其上巷、下巷等不同位置每隔15m或20m設置一個測點，一個工作面會設置十幾至幾十個測點，來測量不同位置的應力情況。可以收集每天整個工作面的電磁輻射最大強度、平均強度和脈沖數3個時間序列，來反映其應力變化。

本文電磁輻射數據來源于躍進煤礦23070工作面。該工作面在2011年2月到2012年1月進行上、巷掘進，2012年2月至12月進行切眼掘進，在切眼掘進期間共發生13次沖擊地壓，而上下巷掘進期間基本未發生沖擊事故。由于不同時期煤巖系統的動力學特性發生了變化，因此將其數據分成兩組，第一組為2011／02／25至2012／01／31上下巷掘進期間的電磁輻射數據，第二組為2012／02／01至2012／12／31切眼掘進期間的數據。

3.3 結合多種監測方式獲得的監測數據

躍進煤礦25110工作面采深達1040m，所采煤層具有強沖擊傾向性，直接頂為厚30m塊狀易碎的泥巖，老頂為礫巖，采用綜放開采方式易導致容易導致應力集中。該工作面沖擊事故多發，僅回采期間就發生10余次沖擊破壞事故。躍進煤礦安裝有ARAMIS微震監測系統，該工作面微震數據完整。另外采用KBD5電磁輻射儀監測獲得了該工作面的電磁輻射數據［5］。作者收集了 2011／07／01 至2012／10／31該工作面回采期間的微震累計能量和頻次、電磁輻射最大強度、平均強度和脈沖數5個時間序列數據。

躍進煤礦23130工作面頂板覆400余米巨厚礫巖，也是一個沖擊地壓事故多發的工作面。本文收集了該工作面2008／10／05至2009／03／20的每天電磁輻射最大強度、平均強度和脈沖數，和0～5m平均煤粉量、6～10m平均煤粉量5個時間序列。

4 沖擊地壓數據的d2計算及結果分析

4.1 微震時間序列的d2計算及結果分析

首先求解平煤十一礦和千秋煤礦微震時間序列的關聯維數。對原始數據歸一化后，采用微震累計能量、頻次2個變量進行相空間重構，重構參數選擇m1＝3，m2＝2，τ1＝1，τ2＝2，然后依次增加m1和m2使m＝5：2：19增加，通過改進的G-P算法得到的十一礦數據的雙對數lnCN（r）-lnr見圖2（a）。同樣的，得到千秋煤礦微震數據的lnCN（r）-lnr曲線，見圖2（b）。

圖2 微震數據的關聯維數求解

在m＝9時，取圖2（a）的線性段lnr＝－0.425～－0.153和圖2（b）的線性段lnr＝－0.362 ～－0.113中的20個點進行最小二乘擬合，得到兩礦微震數據的關聯維數d2分別為2.48和3.92。也就是說，至少需要3個變量才能描述平煤十一礦的微震監控信息，4個變量才能描述躍進礦的微震監控信息，躍進礦的監控信息更加復雜。

在多變量微震時間序列相空間重構的基礎上，利用改進的G-P算法計算得到的躍進煤礦及其25110工作面的關聯維數見表2。可以看出，多變量微震時間序列的關聯維數d2均是分數維數，定性說明微震時間序列具有混沌特性。

比較表2數據，發現平煤十一礦的關聯維數明顯小于千秋和躍進煤礦的結果，表明千秋和躍進煤礦的監測數據復雜程度要強于平煤十一礦。實際上，在監測期間，平煤十一礦發生沖擊事故僅十余次，且破壞性都不大；千秋煤礦和躍進煤礦發生的破壞性沖擊地壓達20多次，其中多次事故造成嚴重的損失和人員傷亡。對比躍進煤礦及其25110工作面的關聯維數發現，全礦微震監測數據的關聯維數明顯大于其下屬工作面的結果，這說明全礦范圍內的沖擊地壓孕育發生過程的復雜程度明顯強于其子工作面。因此，可以得出結論，破壞性沖擊事故越多，破壞程度越強，其關聯維數越大。

表2 多變量微震時間序列的關聯維數

4.2 電磁輻射數據的d2計算及結果分析

23070的工作面的兩組電磁輻射數據的關聯維數時，采用電磁輻射最大強度、平均強度和脈沖數三個變量進行相空間重構后，改進的G-P算法求得兩組數據的關聯維數d2分別為2.43和4.03，都是分數維數，這說明電磁輻射時序也具有混沌特性。對同一個工作面，由于沖擊情況不同，我們看到計算出的關聯維數有明顯差別，上下巷掘進期間，基本沒有沖擊危險，第一組數據的復雜程度不高，關聯維數較小。而在切眼掘進期間，發生了10余次沖擊地壓破壞事故，明顯第二組數據的更加復雜，對應的關聯維數較大。23070工作面兩組數據的結果再次說明，沖擊地壓破壞事故次數越多，破壞程度越強，其關聯維數越大。

4.3 結合多種監測數據的d2計算及結果分析

25110工作面的監測數據包括微震和電磁輻射數據，采用2個微震監測變量重構可以計算微震時序的關聯維數；采用3個電磁輻射時序重構后，得到輻射時序的關聯維數；采用微震和電磁輻射5個時序重構后，得到所有監測時序的關聯維數。結果見表3。類似的，23130工作面的監測時序包括電磁輻射和轉屑量兩類數據，僅采用電磁輻射或鉆屑量一類數據和結合兩類數據計算的關聯維數也見表3。

表3 多類監測數據的關聯維數

比較表3結果我們發現，對同一工作面，不同類型數據計算得到的關聯維數差別不大，都是分數維數。這說明沖擊地壓監測時序具有混沌特性，而關聯維數基本一致的原因是這些監測數據都是該工作面沖擊地壓孕育發展和產生過程中各種信息的反映，同一工作面其過程相同，因此數據復雜程度接近，對應的關聯維數也接近。

5 總結

基于多變量重構理論，對求解單變量時序關聯維數的G-P算法進行擴展和改進以求解多變量時序的關聯維數。以Lorenz系統為例，發現改進G－P算法求解多變量時序的關聯維數時，對數據長度的要求降低、魯棒性強、求解結果更加準確。求得沖擊地壓監測時序的關聯維數均為分數維，說明其具有混沌特性。另外，不同監測數據的關聯維數表明，沖擊地壓事故破壞性越強，其關聯維數越大，監測數據的復雜程度越高。同一工作面的不同種類監測數據的關聯維數接近，說明反映同一工作面沖擊地壓孕育發展、產生過程中各類數據復雜程度接近。這些結論說明可以利用計算的關聯維數值間接判斷沖擊危險性，為沖擊危險預測提供了一種新方法。另外，由于證明監測數據具有混沌特性，這為利用混沌理論預測沖擊危險性提供了依據。

［1］ Xie H P，Chen Z D.Fractal geometry and fracture of rock［J］.Acta Mechanica Sinica，1988，4（3）：255-264.

［2］ 黃滾，尹光志.沖擊地壓粘滑失穩的混沌特性［J］.重慶大學學報，2009，2（6）：633-637.

［3］ 夏永學，潘俊鋒，王元杰，等.基于高精度微震監測的煤巖破裂與應力分布特征研究［J］.煤炭學報，2011，36（2）：239-243.

［4］ 陳國祥，竇林名，曹安業，等.電磁輻射法評定沖擊礦壓危險等級及應用［J］.煤炭學報，2008，33（8）：866-870.

［5］ 陸振裕，竇林名，徐學鋒，等.鉆屑法探測巷道圍巖應力及預測沖擊危險新探究［J］.煤炭工程，2011（1）：71-74.

［6］ 劉貞堂，趙恩來，王恩元，等.不同尺度電磁輻射時間序列的混沌特征初步分析［J］.煤炭學報，2009，34（2）：224-227.

［7］ 陶慧，馬小平，喬美英.基于微震時間序列的沖擊地壓混沌特性分析［J］.煤炭安全，2012，3：52-55.

［8］ Takens F.On the numerical determination of the dimension of an attractor［J］.Lecture notes in Mathematics，1981（898）：361-381.

［9］ 趙善坤，李宏艷，劉軍，等.深部沖擊危險礦井多參量預測預報及解危技術研究［J］.煤炭學報，2011，36（S2）：339-345.

［10］ 王海燕，盧山.非線性時間序列分析及其應用［M］.北京：科學出版社，2006.

［11］ Grassberger P，Procaccia I.Characterization of strange attractors［J］.Physical review letters，1983，50（5）：346-49.

［12］ 王海燕，盛昭瀚，張進.多變量時間序列復雜系統的相空間重構［J］.東南大學學報：自然科學版，2003，33（1）：115-118.

［13］ Lai Y C，Lemer D.Efective scaling regime for computing the correlation dimension from chaotic time series［J］.Physica D，1998（115）：1-18.

Rock burst complexity analysis based on multivariate monitoring time series

TAO Hui1，2，LI Ying3，MA Xiao-ping2
（1.School of Electrical Engineering and Automation，Henan Polytechnic University，Jiaozuo 454000，China；
2.School of Information and Electrical Engineering，China University of Mining and Technology，Xuzhou 221116，China；
3.School of Information Engineering，Jiaozuo University，Jiaozuo 454000，China）

This paper studies the complexity of the data to monitor Rock burst through computing correlation dimension based on phase-space reconstruction of multivariate time series.Given that Rock burst monitoring data had limited-length and contained noise，traditional G-P algorithm is extended and improved.The principle of improved G-P algorithm of solving correlation dimension of multivariate time series was provided，and algorithm was verified through employing it to Lorenz chaotic system.Then a mass of timeseries data to monitor Rock burst were collected by diverse equipment under different burst degree，and their correlation dimensions were computed through improved G-P algorithm.The results demonstrate that the data have chaotic characteristic，and the larger correlation dimension is，the more complex monitoring data is，the stronger Rock-burst damage of corresponding coal mine is.Our achievement can give a novelty approach and basis to predict Rock burst risk based on chaos method.

rock burst；multivariate time-series；correlation dimension；chaotic characteristic

TD324，O415.5

A

1004-4051（2015）10-0118-05

2014-04-23

國家自然科學基金項目“基于多傳感器信息融合技術的礦井瓦斯突出預警研究”資助（編號：60974126）

陶慧（1979－），女，漢族，河南信陽人，副教授，博士，主要研究領域為混沌預測和智能計算。E-mail：65724138＠qq.com。