吳啟紅,萬世明,彭文祥
(1. 成都大學(xué) 城鄉(xiāng)建設(shè)學(xué)院,四川 成都,610106;2. 中南大學(xué) 地球科學(xué)與信息物理學(xué)院,湖南 長沙,410083)
礦產(chǎn)資源在長期的采掘生產(chǎn)過程中,許多礦山形成了數(shù)以百計(jì)的采空區(qū)。特別是一些國有礦山因開采條件復(fù)雜且民采破壞嚴(yán)重, 民間掠奪式采礦留下的群采空區(qū),更是大小不一,形態(tài)各異,層位復(fù)雜,造成上部圍巖多次垮塌,屬于重大事故隱患區(qū)。傳統(tǒng)評價(jià)地下采空區(qū)頂板穩(wěn)定性的方法基本采用半定量分析,包括頂板厚跨比法[1]、魯佩涅依特理論、荷載傳遞線交匯法[2]等。該類方法對頂板的受力情況與破壞機(jī)制認(rèn)識不足,不能全面反映采空區(qū)的應(yīng)力、應(yīng)變分布及破壞狀況,故其計(jì)算結(jié)果的可靠性受到了很大影響,應(yīng)用范圍受到一定限制,為此,許多研究者對其進(jìn)行了研究。如:蔣衛(wèi)東等[3]采用空區(qū)數(shù)量、空區(qū)上部巖體厚度、裂縫密度及裂縫集中度4個(gè)指標(biāo)對采空區(qū)上部地表穩(wěn)定性進(jìn)行灰色定權(quán)聚類分析;李曉靜等[4]綜合考慮圍巖強(qiáng)度和剛度、洞室高度和埋深、水平向地應(yīng)力及其相互間距等影響因素,采用 FLAC3D軟件建立影響地下洞室穩(wěn)定性的層次分析法模型并進(jìn)行分析;凌標(biāo)燦等[5]運(yùn)用人工神經(jīng)網(wǎng)絡(luò)技術(shù),綜合巖石單軸抗壓強(qiáng)度、巖石質(zhì)量指標(biāo)、煤體強(qiáng)度、地下水狀況、工作面月推進(jìn)速度等指標(biāo),建立了采場頂板穩(wěn)定性動(dòng)態(tài)預(yù)測模型用于預(yù)測新集井田頂板穩(wěn)定性分區(qū),并提出不同頂板類型的控制對策;周科平等[6]采用數(shù)值模擬方法分析研究頂板安全厚度與各影響因素之間的關(guān)系,并建立安全頂板厚度非線性神經(jīng)網(wǎng)絡(luò)預(yù)測模型;甄云軍等[7]利用強(qiáng)度折減技術(shù)和二分法原理,以塑性區(qū)的貫通作為頂板破壞的標(biāo)準(zhǔn),得出當(dāng)采空區(qū)頂板的安全系數(shù)為1.2時(shí),各種跨度空區(qū)在不同巖層中的最小安全頂板厚度;彭欣等[8]應(yīng)用 ANSYS三維有限元程序研究了采空區(qū)處理前的圍巖應(yīng)力應(yīng)變情況,分析了采空區(qū)在充填處理后其對近區(qū)開采的影響;何忠明等[9]利用FLAC3D軟件建立雙層空區(qū)數(shù)值計(jì)算模型,根據(jù)厚度折減理論得到安全頂板厚度與空區(qū)跨度之間符合線性關(guān)系,并利用 FLAC3D同時(shí)考慮巖土體的應(yīng)變軟化特性,計(jì)算開挖過程中的變形和破壞情況,分析某大型地下采場開挖穩(wěn)定性;董志宏等[10]針對地下洞室開挖計(jì)算與分析的需要,對已開發(fā)的數(shù)值流形元方法程序進(jìn)行了適當(dāng)改進(jìn)。利用該軟件就某水利樞紐地下洞室穩(wěn)定性問題進(jìn)行了計(jì)算分析,對錨固支護(hù)前后的位移進(jìn)行了比較;李治國[11]利用有限元程序分析了四川省達(dá)縣巴彭公路鐵山下采空區(qū)頂板巖層的變形和破壞機(jī)理,并研究了采空區(qū)治理技術(shù);林杭等[12]借鑒強(qiáng)度折減法計(jì)算邊坡安全系數(shù)的思路,提出采空區(qū)安全頂板預(yù)測的厚度折減法,不斷調(diào)整頂板厚度,直到其達(dá)到臨界破壞狀態(tài),此時(shí)對應(yīng)的頂板厚度即為安全頂板厚度,并對采空區(qū)周圍的應(yīng)力和變形進(jìn)行分析。鑒于目前人們對于2層以上的多層的民采空區(qū)群的穩(wěn)定性未進(jìn)行深入研究,在此,本文作者提出一種“數(shù)值模擬(整體分析)+多級模糊評判(單個(gè)分析)”的綜合評價(jià)方法,并通過工程實(shí)例進(jìn)行分析研究。
礦區(qū)位于廣西壯族自治區(qū),地層以中泥盆統(tǒng)為主,主要分布于礦區(qū)的中部背斜核部,各地層自下而上依次如下。
(1) 堅(jiān)硬的生物礁灰?guī)r巖組(D21)。礦物成分主要為方解石,具生物碎屑結(jié)構(gòu),塊狀構(gòu)造。據(jù)前期礦區(qū)勘查資料,巖石質(zhì)量指標(biāo) RQD為 44%~100%,平均值為75%。根據(jù)井巷水文地質(zhì)工程地質(zhì)調(diào)查,生物礁灰?guī)r巖石質(zhì)量應(yīng)為中等~好,巖體中等完整~較完整。
生物礁灰?guī)r與圍巖常呈指狀交錯(cuò),礁體在構(gòu)造作用下常發(fā)生斷裂破碎,利于礦液的充填,成為本區(qū)的主要容礦空間,形成了大型特富的礦體。
(2) 軟硬相間的泥巖與泥灰?guī)r夾硅質(zhì)巖巖組(D22)、堅(jiān)硬的硅質(zhì)巖巖組(D31)、堅(jiān)硬的扁豆?fàn)罨規(guī)r、條帶狀灰?guī)r巖組(D32)、軟硬相間的泥灰?guī)r、頁巖巖組(D33),環(huán)繞生物礁灰?guī)r巖組(D21)分布,遠(yuǎn)離該礦體,對該礦體開采沒有影響。
(3) 第四系。主要分布于同車江一帶和溝谷中,厚度為3~25 m,賦存了大量砂錫礦。
此外,該礦區(qū)內(nèi)的褶皺及斷裂帶與本次研究的礦體及采空區(qū)距離較遠(yuǎn)。礦區(qū)中部為生物礁灰?guī)r巖溶水區(qū),四周為裂隙水區(qū)。
在礦體內(nèi)井巷及采空區(qū)密集分布。采空區(qū)的疊置關(guān)系見表1,采空區(qū)群的分布見圖1,疊置采空區(qū)見表2。采空區(qū)周邊未發(fā)現(xiàn)有明顯的、新的片落和巖石壓裂現(xiàn)象,沒有巖爆等地壓顯現(xiàn)跡象。
表1 采空區(qū)的疊置關(guān)系Table 1 Overlapping relationship among goafs
圖1 采空區(qū)群的分布圖Fig.1 Distributing graph of goafs
該礦體采空區(qū)群分布復(fù)雜,通過收集礦體采空區(qū)資料,采用 AutoCAD軟件建立三維礦體及采空區(qū)模型,通過“文件-輸出”命令將已建立好的三維礦體及空區(qū)體模型輸出為“*.sat”文件,并導(dǎo)入ANSYS軟件,建立大型復(fù)雜采空區(qū)的計(jì)算模型,最后將ANSYS內(nèi)的結(jié)點(diǎn)和單元信息導(dǎo)入FLAC3D軟件進(jìn)行力學(xué)分析。
礦體復(fù)雜采空區(qū)計(jì)算模型如圖2所示。模型力學(xué)邊界為:在計(jì)算模型左、右和下邊界均為固定位移邊界條件 ux=uy=uz=0,初始應(yīng)力場為巖體的自重應(yīng)力,模型x×y×z=270 m×250 m×490 m,共49萬多個(gè)計(jì)算單元。本次計(jì)算模型中巖性較簡單,僅分為錫礦體和圍巖(灰?guī)r) 2種。計(jì)算所需的巖體物理力學(xué)參數(shù)是根據(jù)室內(nèi)試驗(yàn)確定的巖石物理力學(xué)參數(shù)折減得來的[13],具體參數(shù)結(jié)果見表2。
計(jì)算過程如下:(1) 進(jìn)行自重應(yīng)力場的計(jì)算,生成采空區(qū)的初始應(yīng)力場環(huán)境;(2) 采用摩爾-庫侖模型作為屈服準(zhǔn)則,從上至下分步開挖采空區(qū),將各采空區(qū)設(shè)置為 null模型;(3) 監(jiān)控各采空區(qū)頂板中點(diǎn)的豎向位移;(4) 將FLAC3D軟件分析結(jié)果導(dǎo)入Tecplot進(jìn)行結(jié)果分析,得到各計(jì)算截面上的應(yīng)力分布、塑性區(qū)分布。
表2 采空區(qū)計(jì)算參數(shù)Table 2 Calculation parameters of goafs
通過計(jì)算發(fā)現(xiàn):1號采空區(qū)的頂板及1號和2號采空區(qū)間的重疊頂板按表2進(jìn)行計(jì)算時(shí),整個(gè)空區(qū)群都處于計(jì)算不收斂狀態(tài),這說明1號采空區(qū)的頂板及1號和2號采空區(qū)間的重疊頂板不具備有安全儲備,是極不穩(wěn)定頂板。而事實(shí)上,現(xiàn)場勘察發(fā)現(xiàn)該礦體內(nèi)1號采空區(qū)的頂板及1號和2號采空區(qū)間重疊頂板大部分已經(jīng)塌落。為使計(jì)算能繼續(xù)進(jìn)行,對1號采空區(qū)的頂板及1號和2號采空區(qū)間的重疊頂板的強(qiáng)度逐級增加,直至迭代計(jì)算收斂為止。令1號采空區(qū)的頂板及1號和2號采空區(qū)間的重疊頂板的內(nèi)摩擦角和內(nèi)聚力均取表2中的4.0倍,這樣保證了整個(gè)迭代計(jì)算能繼續(xù)進(jìn)行。
圖2 礦體采空區(qū)群計(jì)算模型及網(wǎng)格劃分Fig.2 Numerical calculation model and mesh generation of orebody and goafs
在計(jì)算模型中,沿Z軸截取若干平面,使各平面通過各采空區(qū)頂板中點(diǎn),得到各采空區(qū)計(jì)算平面的最大、小主應(yīng)力分布(壓應(yīng)力為負(fù),拉應(yīng)力為正)為及塑性區(qū)分布。在各采空區(qū)中,選取相鄰空間區(qū)間頂板厚度小于30 m和頂板跨度較大的空區(qū)群作為研究對象。
2.3.1 多層采空區(qū)計(jì)算平面的應(yīng)力分布及塑性區(qū)分析
沿計(jì)算模型的Z軸截取若干平面,分析計(jì)算平面上各空區(qū)的應(yīng)力分布與塑性區(qū)分布,尤其是多層空區(qū)間頂板的應(yīng)力與塑性區(qū)分布。這里僅以1號、2號、13號和21號采空區(qū)群為例說明分析過程。
圖3 1號、2號、13號和21號空區(qū)相對位置Fig.3 Relative positions of No. 1, 2, 13 and 21
圖4 計(jì)算截面上的最大主應(yīng)力分布(Z=108.03 m)Fig.4 Distributing of maximum principal stress (Z=108.03 m)
2號采空區(qū)面積為3 632 m2,高為2 m,跨度達(dá)30 m左右。2號采空區(qū)頂板距1號采空區(qū)底板僅2 m,中部頂板部分垮下。21號采空區(qū)位于礦體標(biāo)高為-115~-120 m,面積為884 m2,高為5.0 m,與13號采空區(qū)和2號采空區(qū)部分疊置。圖3所示為1號、2號、13號和21號空區(qū)相對位置圖。取Z=108.038 5的截面進(jìn)分析,計(jì)算截面上的最大主應(yīng)力分布(圖4)。從圖4可見:1號、2號、21號和13號空區(qū)間頂板均為卸荷區(qū)且最小主應(yīng)力為拉應(yīng)力,為0.4~0.8 MPa(圖5)。圖6所示為1號、2號、21號和13號采空區(qū)間頂板的塑性區(qū)分布。從圖6可見:1號和2號采空區(qū)間頂板全部處于塑性狀態(tài);1號和2號空區(qū)間的頂板很薄,其厚度僅為2 m左右,該頂板處于極不穩(wěn)定狀態(tài);21號采空區(qū)和2號采空區(qū)間頂板塑性區(qū)范圍大而且接近貫通,21號采空區(qū)頂板在計(jì)算截面上較不穩(wěn)定;而13號采空區(qū),頂板塑性區(qū)范圍較少,且13號采空區(qū)頂板和21號采空區(qū)底板相距11 m,可以認(rèn)為13號空區(qū)在計(jì)算截面上是穩(wěn)定的。
圖5 1號、2號、13號和21號采空區(qū)間頂板的最小主應(yīng)力分布Fig.5 Distribution of minimum principal stress of Nos 1, 2, 13 and 21
圖6 1號、2號、13號和21號采空區(qū)間頂板的塑性區(qū)分布Fig.6 Distribution of plastic zone of No. 1, 2, 13 and 21
2.3.2 各采空頂板中點(diǎn)位移監(jiān)控
在計(jì)算過程中監(jiān)控各采空區(qū)頂板中點(diǎn)的豎向位移,得到各采空區(qū)中點(diǎn)的豎向位移隨計(jì)算時(shí)步的演化曲線。研究結(jié)果表明:2號結(jié)點(diǎn)為2號采空區(qū)頂板中點(diǎn)的豎向位移曲線,盡管在計(jì)算中將1號采空區(qū)及1號和2號采空區(qū)間的重疊頂板的強(qiáng)度增加了4.0倍,但2號采空區(qū)頂板中點(diǎn)的豎向位移不穩(wěn)定發(fā)展,最終位移達(dá)6.46 cm,該頂板不穩(wěn)定;1號采空區(qū)頂板中點(diǎn)的豎向位移達(dá)3.62 cm,這2個(gè)采空區(qū)跨度大,頂板存在過大的塑性變形,因此,1號和2號采空區(qū)頂板極為不穩(wěn)定。其他各采空區(qū)頂板中點(diǎn)的豎向位移為1.88~2.90 cm。這表明:總體上礦體復(fù)雜采空區(qū)頂板存在較大的塑性變形,各采空區(qū)頂板位移比較大;1號和2號采空區(qū)頂板(尤其是2號)出現(xiàn)了大規(guī)模的塑性流動(dòng)。
對所建立的礦體采空區(qū)群模型計(jì)算結(jié)果從應(yīng)力分布、塑性區(qū)以及各采空區(qū)頂板中點(diǎn)位移進(jìn)行分析,各采空區(qū)穩(wěn)定性分類結(jié)果見表3(穩(wěn)定性分類:III為不穩(wěn)定采空區(qū),II為較不穩(wěn)定采空區(qū),I為穩(wěn)定采空區(qū))。
礦體采空區(qū)穩(wěn)定性模糊評價(jià)采用的是吳啟紅等[14]前期所建立的采空區(qū)穩(wěn)定性多級模糊評價(jià)模型。該多級模糊評判模型可針對各類復(fù)雜的采空區(qū),通過選取不同的因素、因子對其進(jìn)行穩(wěn)定性評價(jià),具有較強(qiáng)的實(shí)用性和靈活性;同時(shí),為了保證各因子具有等效性和同序性,在進(jìn)行模糊運(yùn)算之前,首先對因子原始數(shù)據(jù)進(jìn)行處理,即對因子的各級界限值進(jìn)行標(biāo)準(zhǔn)化處理。并由于預(yù)測因子有定性(離散型變量)和定量(連續(xù)型變量)指標(biāo)2類,對于離散型變量,隸屬度取值采用專家評定法;而對連續(xù)型變量的隸屬度取值,通過代入實(shí)測值采用三相線性隸屬函數(shù)計(jì)算得到。同時(shí),對因子的權(quán)重采用計(jì)算機(jī)反演的方法來確定。
表3 礦體各采空區(qū)情況和穩(wěn)定性分類Table 3 Classification of stability and condition of each goaf
運(yùn)用該模糊評價(jià)模型的評價(jià)對礦體各采空區(qū)的穩(wěn)定性進(jìn)行分類,見表3。
通過對比多層采空區(qū)群三維數(shù)值模擬分析和多級模糊評判的結(jié)果可以發(fā)現(xiàn):兩者對于該礦體的采空區(qū)穩(wěn)定性分類基本一致,僅在13號采空區(qū)的穩(wěn)定性分類上存在差異(數(shù)值模擬巖石力學(xué)分析結(jié)果為I,即穩(wěn)定;多級模糊評判結(jié)果為II,即較不穩(wěn)定)。但經(jīng)綜合分析后評價(jià)其穩(wěn)定性為II,即為較不穩(wěn)定。其原因有以下2點(diǎn):(1) 13號采空區(qū)在三維上為不規(guī)則形狀,而在本次數(shù)值模擬分析時(shí)所取的計(jì)算截面(Z=108.03 m)并不一定為該采空區(qū)最不利的計(jì)算截面;(2) 目前該采空區(qū)或許處于相對穩(wěn)定狀態(tài),在沒有達(dá)到巖體的破壞極限之前,其相對穩(wěn)定狀態(tài)仍然會(huì)持續(xù)下去,當(dāng)受到新的擾動(dòng)時(shí),極可能引起空區(qū)坍塌、冒頂、圍巖開裂或?qū)娱g錯(cuò)動(dòng),甚至帶來危害較大的地壓活動(dòng)。為避免上述危害,從治理方案的選取上考慮,評判其為較不穩(wěn)定。
(1) 通過AUTOCAD和ANSYS建立某礦體復(fù)雜多層采空區(qū)群三維模擬計(jì)算模型,并利用 FLAC3D軟件從采空區(qū)頂板應(yīng)力場、塑性區(qū)及頂板中點(diǎn)位移演化幾方面對其進(jìn)行巖石力學(xué)穩(wěn)定性分析。該方法為類似復(fù)雜多層采空區(qū)群數(shù)值模型的建立以及穩(wěn)定性分析提供了一種簡單易行的思路。
(2) 所建立的三維數(shù)值模擬分析模型,通過選取不同的截面(Z軸不同),采空區(qū)頂板在不同部位顯示出的穩(wěn)定性不一樣,例如19號采空區(qū)頂板在不同計(jì)算截面上的穩(wěn)性相差很大,這一點(diǎn)與現(xiàn)場勘測得到的結(jié)果是一致的。
(3) 通過數(shù)值模擬分析以及多級模糊評判這 2種方法對同一礦體各采空區(qū)群穩(wěn)定性進(jìn)行分類比較發(fā)現(xiàn),僅在13號采空區(qū)的穩(wěn)定性分類上存在較小差異,這也從側(cè)面驗(yàn)證了該2種方法具有較高的準(zhǔn)確度。
(4) 該綜合評價(jià)法通過多種方法從整體到單個(gè)的對比分析提高了多層采空區(qū)群穩(wěn)定性分析的準(zhǔn)確度,這為后續(xù)的礦山安全治理研究提供了更加可靠的依據(jù)。
[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] 蔣衛(wèi)東, 李夕兵, 胡柳青, 等. 基于灰色定權(quán)聚類的采空區(qū)上部地表穩(wěn)定性分析[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] 李曉靜, 朱維申, 陳衛(wèi)忠, 等. 層次分析法確定影響地下洞室圍巖穩(wěn)定性各因素的權(quán)值[J]. 巖石力學(xué)與工程學(xué)報(bào), 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] 凌標(biāo)燦, 彭蘇萍, 張滇河, 等. 采場頂板穩(wěn)定性動(dòng)態(tài)工程分類[J]. 巖石力學(xué)與工程學(xué)報(bào), 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] 周科平, 蘇家紅, 古德生, 等. 復(fù)雜充填體下礦體開采安全頂板厚度非線性預(yù)測方法[J]. 中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 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] 甄云軍, 陳開翔, 劉應(yīng)發(fā), 等. 地下采空區(qū)頂板安全厚度的確定[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] 彭欣, 崔棟梁, 李夕兵, 等. 特大采空區(qū)近區(qū)開采的穩(wěn)定性分析[J]. 中國礦業(yè), 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] 何忠明, 彭振斌, 曹平, 等. 雙層采空區(qū)開挖頂板穩(wěn)定性的FLAC3D數(shù)值分析[J]. 中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 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] 董志宏, 鄔愛清,丁秀麗. 基于數(shù)值流形元方法的地下洞室穩(wěn)定性分析[J]. 巖石力學(xué)與工程學(xué)報(bào), 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] 李治國. 鐵山隧道采空區(qū)穩(wěn)定性分析及治理技術(shù)研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 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] 林杭, 曹平, 李江騰, 等. 采空區(qū)臨界安全頂板預(yù)測的厚度折減法[J]. 煤炭學(xué)報(bào), 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] 《工程地質(zhì)手冊》編寫委員會(huì). 工程地質(zhì)手冊[M]. 4版. 北京:中國建筑工業(yè)出版社, 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] 吳啟紅, 彭振斌, 陳科平, 等. 礦山采空區(qū)穩(wěn)定性二級模糊綜合評判. 中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 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.