一種多層采空區群穩定性的綜合評價法

吳啟紅，萬世明，彭文祥

(1. 成都大學 城鄉建設學院，四川 成都，610106；2. 中南大學 地球科學與信息物理學院，湖南 長沙，410083)

礦產資源在長期的采掘生產過程中，許多礦山形成了數以百計的采空區。特別是一些國有礦山因開采條件復雜且民采破壞嚴重， 民間掠奪式采礦留下的群采空區，更是大小不一，形態各異，層位復雜，造成上部圍巖多次垮塌，屬于重大事故隱患區。傳統評價地下采空區頂板穩定性的方法基本采用半定量分析，包括頂板厚跨比法[1]、魯佩涅依特理論、荷載傳遞線交匯法[2]等。該類方法對頂板的受力情況與破壞機制認識不足，不能全面反映采空區的應力、應變分布及破壞狀況，故其計算結果的可靠性受到了很大影響，應用范圍受到一定限制，為此，許多研究者對其進行了研究。如：蔣衛東等[3]采用空區數量、空區上部巖體厚度、裂縫密度及裂縫集中度4個指標對采空區上部地表穩定性進行灰色定權聚類分析；李曉靜等[4]綜合考慮圍巖強度和剛度、洞室高度和埋深、水平向地應力及其相互間距等影響因素，采用 FLAC3D軟件建立影響地下洞室穩定性的層次分析法模型并進行分析；凌標燦等[5]運用人工神經網絡技術，綜合巖石單軸抗壓強度、巖石質量指標、煤體強度、地下水狀況、工作面月推進速度等指標，建立了采場頂板穩定性動態預測模型用于預測新集井田頂板穩定性分區，并提出不同頂板類型的控制對策；周科平等[6]采用數值模擬方法分析研究頂板安全厚度與各影響因素之間的關系，并建立安全頂板厚度非線性神經網絡預測模型；甄云軍等[7]利用強度折減技術和二分法原理，以塑性區的貫通作為頂板破壞的標準，得出當采空區頂板的安全系數為1.2時，各種跨度空區在不同巖層中的最小安全頂板厚度；彭欣等[8]應用 ANSYS三維有限元程序研究了采空區處理前的圍巖應力應變情況，分析了采空區在充填處理后其對近區開采的影響；何忠明等[9]利用FLAC3D軟件建立雙層空區數值計算模型，根據厚度折減理論得到安全頂板厚度與空區跨度之間符合線性關系，并利用 FLAC3D同時考慮巖土體的應變軟化特性，計算開挖過程中的變形和破壞情況，分析某大型地下采場開挖穩定性；董志宏等[10]針對地下洞室開挖計算與分析的需要，對已開發的數值流形元方法程序進行了適當改進。利用該軟件就某水利樞紐地下洞室穩定性問題進行了計算分析，對錨固支護前后的位移進行了比較；李治國[11]利用有限元程序分析了四川省達縣巴彭公路鐵山下采空區頂板巖層的變形和破壞機理，并研究了采空區治理技術；林杭等[12]借鑒強度折減法計算邊坡安全系數的思路，提出采空區安全頂板預測的厚度折減法，不斷調整頂板厚度，直到其達到臨界破壞狀態，此時對應的頂板厚度即為安全頂板厚度，并對采空區周圍的應力和變形進行分析。鑒于目前人們對于2層以上的多層的民采空區群的穩定性未進行深入研究，在此，本文作者提出一種“數值模擬(整體分析)+多級模糊評判(單個分析)”的綜合評價方法，并通過工程實例進行分析研究。

1 工程概況

1.1 工程地質條件

礦區位于廣西壯族自治區，地層以中泥盆統為主，主要分布于礦區的中部背斜核部，各地層自下而上依次如下。

(1) 堅硬的生物礁灰巖巖組(D21)。礦物成分主要為方解石，具生物碎屑結構，塊狀構造。據前期礦區勘查資料，巖石質量指標 RQD為 44%～100%，平均值為75%。根據井巷水文地質工程地質調查，生物礁灰巖巖石質量應為中等～好，巖體中等完整～較完整。

生物礁灰巖與圍巖常呈指狀交錯，礁體在構造作用下常發生斷裂破碎，利于礦液的充填，成為本區的主要容礦空間，形成了大型特富的礦體。

(2) 軟硬相間的泥巖與泥灰巖夾硅質巖巖組(D22)、堅硬的硅質巖巖組(D31)、堅硬的扁豆狀灰巖、條帶狀灰巖巖組(D32)、軟硬相間的泥灰巖、頁巖巖組(D33)，環繞生物礁灰巖巖組(D21)分布，遠離該礦體，對該礦體開采沒有影響。

(3) 第四系。主要分布于同車江一帶和溝谷中，厚度為3～25 m，賦存了大量砂錫礦。

此外，該礦區內的褶皺及斷裂帶與本次研究的礦體及采空區距離較遠。礦區中部為生物礁灰巖巖溶水區，四周為裂隙水區。

1.2 現狀采空區圍巖穩定性評價

在礦體內井巷及采空區密集分布。采空區的疊置關系見表1，采空區群的分布見圖1，疊置采空區見表2。采空區周邊未發現有明顯的、新的片落和巖石壓裂現象，沒有巖爆等地壓顯現跡象。

表1 采空區的疊置關系Table 1 Overlapping relationship among goafs

圖1 采空區群的分布圖Fig.1 Distributing graph of goafs

2 礦體多層采空區群穩定性的數值模擬分析

2.1 建模過程

該礦體采空區群分布復雜，通過收集礦體采空區資料，采用 AutoCAD軟件建立三維礦體及采空區模型，通過“文件－輸出”命令將已建立好的三維礦體及空區體模型輸出為“*.sat”文件，并導入ANSYS軟件，建立大型復雜采空區的計算模型，最后將ANSYS內的結點和單元信息導入FLAC3D軟件進行力學分析。

2.2 礦體內采空區群數值分析模型

礦體復雜采空區計算模型如圖2所示。模型力學邊界為：在計算模型左、右和下邊界均為固定位移邊界條件 ux=uy=uz=0，初始應力場為巖體的自重應力，模型x×y×z=270 m×250 m×490 m，共49萬多個計算單元。本次計算模型中巖性較簡單，僅分為錫礦體和圍巖(灰巖) 2種。計算所需的巖體物理力學參數是根據室內試驗確定的巖石物理力學參數折減得來的[13]，具體參數結果見表2。

計算過程如下：(1) 進行自重應力場的計算，生成采空區的初始應力場環境；(2) 采用摩爾-庫侖模型作為屈服準則，從上至下分步開挖采空區，將各采空區設置為 null模型；(3) 監控各采空區頂板中點的豎向位移；(4) 將FLAC3D軟件分析結果導入Tecplot進行結果分析，得到各計算截面上的應力分布、塑性區分布。

表2 采空區計算參數Table 2 Calculation parameters of goafs

通過計算發現：1號采空區的頂板及1號和2號采空區間的重疊頂板按表2進行計算時，整個空區群都處于計算不收斂狀態，這說明1號采空區的頂板及1號和2號采空區間的重疊頂板不具備有安全儲備，是極不穩定頂板。而事實上，現場勘察發現該礦體內1號采空區的頂板及1號和2號采空區間重疊頂板大部分已經塌落。為使計算能繼續進行，對1號采空區的頂板及1號和2號采空區間的重疊頂板的強度逐級增加，直至迭代計算收斂為止。令1號采空區的頂板及1號和2號采空區間的重疊頂板的內摩擦角和內聚力均取表2中的4.0倍，這樣保證了整個迭代計算能繼續進行。

圖2 礦體采空區群計算模型及網格劃分Fig.2 Numerical calculation model and mesh generation of orebody and goafs

2.3 數值模擬結果分析

在計算模型中，沿Z軸截取若干平面，使各平面通過各采空區頂板中點，得到各采空區計算平面的最大、小主應力分布(壓應力為負，拉應力為正)為及塑性區分布。在各采空區中，選取相鄰空間區間頂板厚度小于30 m和頂板跨度較大的空區群作為研究對象。

2.3.1 多層采空區計算平面的應力分布及塑性區分析

沿計算模型的Z軸截取若干平面，分析計算平面上各空區的應力分布與塑性區分布，尤其是多層空區間頂板的應力與塑性區分布。這里僅以1號、2號、13號和21號采空區群為例說明分析過程。

圖3 1號、2號、13號和21號空區相對位置Fig.3 Relative positions of No. 1, 2, 13 and 21

圖4 計算截面上的最大主應力分布(Z=108.03 m)Fig.4 Distributing of maximum principal stress (Z=108.03 m)

2號采空區面積為3 632 m2，高為2 m，跨度達30 m左右。2號采空區頂板距1號采空區底板僅2 m，中部頂板部分垮下。21號采空區位于礦體標高為-115～-120 m，面積為884 m2，高為5.0 m，與13號采空區和2號采空區部分疊置。圖3所示為1號、2號、13號和21號空區相對位置圖。取Z=108.038 5的截面進分析，計算截面上的最大主應力分布(圖4)。從圖4可見：1號、2號、21號和13號空區間頂板均為卸荷區且最小主應力為拉應力，為0.4～0.8 MPa(圖5)。圖6所示為1號、2號、21號和13號采空區間頂板的塑性區分布。從圖6可見：1號和2號采空區間頂板全部處于塑性狀態；1號和2號空區間的頂板很薄，其厚度僅為2 m左右，該頂板處于極不穩定狀態；21號采空區和2號采空區間頂板塑性區范圍大而且接近貫通，21號采空區頂板在計算截面上較不穩定；而13號采空區，頂板塑性區范圍較少，且13號采空區頂板和21號采空區底板相距11 m，可以認為13號空區在計算截面上是穩定的。

圖5 1號、2號、13號和21號采空區間頂板的最小主應力分布Fig.5 Distribution of minimum principal stress of Nos 1, 2, 13 and 21

圖6 1號、2號、13號和21號采空區間頂板的塑性區分布Fig.6 Distribution of plastic zone of No. 1, 2, 13 and 21

2.3.2 各采空頂板中點位移監控

在計算過程中監控各采空區頂板中點的豎向位移，得到各采空區中點的豎向位移隨計算時步的演化曲線。研究結果表明：2號結點為2號采空區頂板中點的豎向位移曲線，盡管在計算中將1號采空區及1號和2號采空區間的重疊頂板的強度增加了4.0倍，但2號采空區頂板中點的豎向位移不穩定發展，最終位移達6.46 cm，該頂板不穩定；1號采空區頂板中點的豎向位移達3.62 cm，這2個采空區跨度大，頂板存在過大的塑性變形，因此，1號和2號采空區頂板極為不穩定。其他各采空區頂板中點的豎向位移為1.88～2.90 cm。這表明：總體上礦體復雜采空區頂板存在較大的塑性變形，各采空區頂板位移比較大；1號和2號采空區頂板(尤其是2號)出現了大規模的塑性流動。

2.4 礦體采空區穩定性分類

對所建立的礦體采空區群模型計算結果從應力分布、塑性區以及各采空區頂板中點位移進行分析，各采空區穩定性分類結果見表3(穩定性分類：III為不穩定采空區，II為較不穩定采空區，I為穩定采空區)。

3 礦體各采空區穩定性模糊評價

礦體采空區穩定性模糊評價采用的是吳啟紅等[14]前期所建立的采空區穩定性多級模糊評價模型。該多級模糊評判模型可針對各類復雜的采空區，通過選取不同的因素、因子對其進行穩定性評價，具有較強的實用性和靈活性；同時，為了保證各因子具有等效性和同序性，在進行模糊運算之前，首先對因子原始數據進行處理，即對因子的各級界限值進行標準化處理。并由于預測因子有定性(離散型變量)和定量(連續型變量)指標2類，對于離散型變量，隸屬度取值采用專家評定法；而對連續型變量的隸屬度取值，通過代入實測值采用三相線性隸屬函數計算得到。同時，對因子的權重采用計算機反演的方法來確定。

表3 礦體各采空區情況和穩定性分類Table 3 Classification of stability and condition of each goaf

運用該模糊評價模型的評價對礦體各采空區的穩定性進行分類，見表3。

4 對比分析

通過對比多層采空區群三維數值模擬分析和多級模糊評判的結果可以發現：兩者對于該礦體的采空區穩定性分類基本一致，僅在13號采空區的穩定性分類上存在差異(數值模擬巖石力學分析結果為I，即穩定；多級模糊評判結果為II，即較不穩定)。但經綜合分析后評價其穩定性為II，即為較不穩定。其原因有以下2點：(1) 13號采空區在三維上為不規則形狀，而在本次數值模擬分析時所取的計算截面(Z=108.03 m)并不一定為該采空區最不利的計算截面；(2) 目前該采空區或許處于相對穩定狀態，在沒有達到巖體的破壞極限之前，其相對穩定狀態仍然會持續下去，當受到新的擾動時，極可能引起空區坍塌、冒頂、圍巖開裂或層間錯動，甚至帶來危害較大的地壓活動。為避免上述危害，從治理方案的選取上考慮，評判其為較不穩定。

5 結論

(1) 通過AUTOCAD和ANSYS建立某礦體復雜多層采空區群三維模擬計算模型，并利用 FLAC3D軟件從采空區頂板應力場、塑性區及頂板中點位移演化幾方面對其進行巖石力學穩定性分析。該方法為類似復雜多層采空區群數值模型的建立以及穩定性分析提供了一種簡單易行的思路。

(2) 所建立的三維數值模擬分析模型，通過選取不同的截面(Z軸不同)，采空區頂板在不同部位顯示出的穩定性不一樣，例如19號采空區頂板在不同計算截面上的穩性相差很大，這一點與現場勘測得到的結果是一致的。

(3) 通過數值模擬分析以及多級模糊評判這 2種方法對同一礦體各采空區群穩定性進行分類比較發現，僅在13號采空區的穩定性分類上存在較小差異，這也從側面驗證了該2種方法具有較高的準確度。

(4) 該綜合評價法通過多種方法從整體到單個的對比分析提高了多層采空區群穩定性分析的準確度，這為后續的礦山安全治理研究提供了更加可靠的依據。

[1] Swift G M, Reddish D J. Stability problems associated with an abandoned ironstone mine[J]. Bulletin of Engineering Geology and the Environment, 2002, 61(3): 227-239.

[2] Nomikos P P, Sofianos A I, Tsoutrelis C E. Structural response of vertically multi-jointed roof rock beams[J]. International Journal of Rock Mechanics and Mining Sciences, 2002, 39(1): 79-94.

[3] 蔣衛東, 李夕兵, 胡柳青, 等. 基于灰色定權聚類的采空區上部地表穩定性分析[J]. 礦冶工程, 2002, 22(4): 15-17.JIANG Wei-dong, LI Xi-bing, HU Liu-qing, et al. Stability analysis of ground above excavated area based on grey fixed weight cluster[J]. Mining and Metallurgical Engineering, 2002,22(4): 15-17.

[4] 李曉靜, 朱維申, 陳衛忠, 等. 層次分析法確定影響地下洞室圍巖穩定性各因素的權值[J]. 巖石力學與工程學報, 2004,23(S2): 4731-4734.LI Xiao-jing, ZHU Wei-shen, CHEN Wei-zhong, et al.Determining weight of factors in stability analysis of underground caverns by analytic hierarchy process[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(S2):4731-4734.

[5] 凌標燦, 彭蘇萍, 張滇河, 等. 采場頂板穩定性動態工程分類[J]. 巖石力學與工程學報, 2003, 22(9): 1474-1477.LING Biao-can, PENG Su-ping, ZHANG Sheng-he, et al.Dynamic engineering classification of stope roof stability[J].Chinese Journal of Rock Mechanics and Engineering, 2003,22(9): 1474-1477.

[6] 周科平, 蘇家紅, 古德生, 等. 復雜充填體下礦體開采安全頂板厚度非線性預測方法[J]. 中南大學學報: 自然科學版, 2005,36(6): 1094-1099.ZHOU Ke-ping, SU Jia-hong, GU De-sheng, et al. The nonlinear forecasting method of the least security coping thickness when mining under complex filling body[J]. Journal of Central South University: Science and Technology, 2005, 36(6): 1094-1099.

[7] 甄云軍, 陳開翔, 劉應發, 等. 地下采空區頂板安全厚度的確定[J]. 化工礦物與加工, 2007(9): 19-20.ZHEN Yun-jun, CHEN Kai-xiang, LIU Ying-fa, et al.Determination of roof safety thickness for underground mined-out area[J]. Industrial Minerals & Processing, 2007(9):19-20

[8] 彭欣, 崔棟梁, 李夕兵, 等. 特大采空區近區開采的穩定性分析[J]. 中國礦業, 2007, 16(4): 70-73.PENG Xin, CUI Dong-liang, LI Xi-bing, et al. Adopt the stability analysis that the empty area near area mines greatly especially[J]. China Mining Magazine, 2007, 16(4): 70-73.

[9] 何忠明, 彭振斌, 曹平, 等. 雙層采空區開挖頂板穩定性的FLAC3D數值分析[J]. 中南大學學報: 自然科學版, 2009,40(4): 1066-1071.HE Zhong-ming, PENG Zhen-bin, CAO Ping, et al. Numerical analysis for roof stability of double gob area after excavation by FLAC3D[J]. Journal of Central South University: Science and Technology, 2009, 40(4): 1066-1071.

[10] 董志宏, 鄔愛清,丁秀麗. 基于數值流形元方法的地下洞室穩定性分析[J]. 巖石力學與工程學報, 2004, 23(S2): 4956-4959.DONG Zhi-hong, WU Ai-qing, DING Xiu-li. Stability study of underground cavern with numerical manifold method[J].Chinese Journal of Rock Mechanics and Engineering, 2004,23(S2): 4956-4959.

[11] 李治國. 鐵山隧道采空區穩定性分析及治理技術研究[J]. 巖石力學與工程學報, 2002, 21(8): 1168-1173.LI Zhi-guo. Stability analysis and reinforcement technology of mined-out area in tieshan tunnel[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(8): 1168-1173.

[12] 林杭, 曹平, 李江騰, 等. 采空區臨界安全頂板預測的厚度折減法[J]. 煤炭學報, 2009, 34(1): 53-57.LIN Hang, CAO Ping, LI Jiang-teng, et al. The thickness reduction method in forecasting the critical safety roof thickness of gob area[J]. Journal of China Coal Society, 2009, 34(1):53-57.

[13] 《工程地質手冊》編寫委員會. 工程地質手冊[M]. 4版. 北京:中國建筑工業出版社, 2007: 165-170.The Compilation Committee of Engineering Geology Handbook.Engineering geology handbook[M]. 4th ed. Beijing: China Architecture and Building Press, 2007: 165-170.

[14] 吳啟紅, 彭振斌, 陳科平, 等. 礦山采空區穩定性二級模糊綜合評判. 中南大學學報: 自然科學版, 2010, 41(2): 661-667.WU Qi-hong, PENG Zhen-bin, CHEN Ke-ping, et al. Synthetic judgment on two-stage fuzzy of stability of mine gob area[J].Journal of Central South University: Science and Technology,2010, 41(2): 661-667.