基于加權(quán)信息法的大孤山露天礦滑坡風(fēng)險(xiǎn)分析

李 旭 于慶磊 楊天鴻 董 鑫 鄧文學(xué)

（東北大學(xué)資源與土木工程學(xué)院，遼寧 沈陽(yáng) 110819）

大孤山露天礦目前已進(jìn)入深部開(kāi)采階段，由于開(kāi)采期間礦山巖體處于動(dòng)態(tài)演化中，導(dǎo)致邊坡的穩(wěn)定性受到諸多因素的影響。為此，有必要進(jìn)行露天礦山的邊坡穩(wěn)定性研究?；谒惴A(yù)測(cè)邊坡發(fā)生失穩(wěn)的高風(fēng)險(xiǎn)地區(qū)，對(duì)高風(fēng)險(xiǎn)進(jìn)行重點(diǎn)加固防護(hù)，這樣有利于提高礦山生產(chǎn)的安全性和節(jié)省邊坡加固的成本。

近年來(lái)，在預(yù)測(cè)滑坡風(fēng)險(xiǎn)的案例中，許多學(xué)者在預(yù)測(cè)滑坡失穩(wěn)方面做出重要貢獻(xiàn)。樊曉一等［1］基于三峽重點(diǎn)滑坡區(qū)域利用層次分析法和地形判別法，建立典型滑坡危險(xiǎn)度判別矩陣進(jìn)行危險(xiǎn)度計(jì)算；白宇等［2］在概率論與數(shù)理統(tǒng)計(jì)方面與預(yù)測(cè)邊坡的極限平衡原理相結(jié)合，利用Monte-Carlo法來(lái)模擬隨機(jī)變量分布；李雪平等［3］基于歷史滑坡案例建立了Logistic回歸模型，利用ArcGIS進(jìn)行影響因子提取，基于VisualC++6.0進(jìn)行Logistic回歸模型分析；王萌等［4］以沐川縣的滑坡為對(duì)象，提取5個(gè)滑坡因子，采用貢獻(xiàn)率法計(jì)算權(quán)重；畢華興等［5］利用航片和TM數(shù)據(jù)提取并構(gòu)建影響滑坡的8個(gè)因子，并建立滑坡因子專題圖，以8張滑坡因子滑坡專題圖為基礎(chǔ)進(jìn)行抽樣調(diào)查，利用數(shù)量化理論公式進(jìn)行數(shù)據(jù)處理，根據(jù)樣本得到的分值進(jìn)行滑坡預(yù)測(cè)；于進(jìn)慶等［6］基于GIS系統(tǒng)提取影響滑坡因子，并建立層次分析法的判斷矩陣，利用GIS統(tǒng)計(jì)獲取的滑坡因子專題圖面積得到的信息量作為層次分析法判別的重要程度進(jìn)行滑坡因子計(jì)算；倪恒等［7］提出用敏感性分析思想來(lái)計(jì)算各因素對(duì)滑坡產(chǎn)生影響的大小，以趙樹(shù)嶺滑坡為實(shí)驗(yàn)對(duì)象，考慮了5個(gè)影響因素，通過(guò)滑坡敏感性分析結(jié)合SARMA算法得出結(jié)論：水位變化＞地震荷載＞內(nèi)摩擦角＞粘聚力＞巖體密度；張桂榮等［8］依托GIS平臺(tái)，采用規(guī)則網(wǎng)格單元?jiǎng)澐址椒?，運(yùn)用信息量模型對(duì)該區(qū)斜坡穩(wěn)定性進(jìn)行了空間定量預(yù)測(cè)，并依信息量法的結(jié)果編制了該區(qū)的危險(xiǎn)性預(yù)測(cè)分區(qū)圖；孫冉、李娜娜、吳麗清等［9-11］學(xué)者基于加權(quán)信息量算法，結(jié)合當(dāng)?shù)氐刭|(zhì)條件提取影響因素，并通過(guò)GIS軟件制作滑坡風(fēng)險(xiǎn)云圖。

綜上分析，眾學(xué)者已對(duì)滑坡風(fēng)險(xiǎn)預(yù)測(cè)做了諸多貢獻(xiàn)，但都是基于區(qū)域大、滑坡點(diǎn)多的情況。對(duì)于在礦山中使用這類預(yù)測(cè)方法并不多見(jiàn)。故在本項(xiàng)目中，將加權(quán)信息量法運(yùn)用在礦山邊坡風(fēng)險(xiǎn)預(yù)測(cè)方面。通過(guò)對(duì)大孤山露天礦的數(shù)字高程模型（DEM）進(jìn)行空間分析和地形數(shù)據(jù)提取，利用GIS的柵格圖繪制和矢量線段構(gòu)建緩沖區(qū)功能提取礦山地質(zhì)巖性和斷層因子?；趯哟畏治龇ńo諸多因素進(jìn)行重要性評(píng)價(jià)，并用GIS柵格統(tǒng)計(jì)進(jìn)行信息量的評(píng)估，最終得出滑坡風(fēng)險(xiǎn)圖，從而達(dá)到提高礦山施工安全性，節(jié)約加固資金成本的目的。

1 大孤山鐵礦地質(zhì)概況

大孤山鐵礦位于鞍山市東南12 km的千山腳下，礦區(qū)中心點(diǎn)坐標(biāo)：東經(jīng) 123°3'35.2443"、北緯41°3'5.22813929688"，占地10.6 km2，是亞洲最深的露天鐵礦。礦區(qū)主要巖性為太古界鞍山群櫻桃園組和元古界遼河群浪子山組的中深變質(zhì)巖系，其次為震旦系釣魚(yú)臺(tái)石英巖及第四紀(jì)地層［12］。目前礦山（圖1）已開(kāi)采至-330 m水平，開(kāi)采深度420 m［13］。由于礦山進(jìn)入深部開(kāi)采階段，滑坡現(xiàn)象及巖石開(kāi)裂變形情況愈加嚴(yán)重，引起大孤山鐵礦破壞的原因［14］主要有：結(jié)構(gòu)面和人工開(kāi)挖爆破?，F(xiàn)階段，大孤山西北端已出現(xiàn)多處裂縫，且仍處在變形階段，嚴(yán)重威脅巷道運(yùn)輸和礦山作業(yè)安全。

根據(jù)現(xiàn)場(chǎng)地質(zhì)構(gòu)造、水文地質(zhì)和巖性等諸多條件，將大孤山礦坑已發(fā)生滑坡地區(qū)分為3個(gè)區(qū)域。其中I區(qū)位于西井邊坡F15斷層南部，巖性為太古代花崗巖（糜棱巖）；II區(qū)位于F15斷層以北、F14斷層以南，受2條斷層帶影響呈楔形發(fā)育；III區(qū)位于F14斷層以北，其巖性也為太古代花崗巖。邊坡破壞現(xiàn)狀如圖2所示。

2 滑坡評(píng)價(jià)因子的選取和層次分析法權(quán)重的確定

2.1 滑坡評(píng)價(jià)因子的選取

考慮當(dāng)前礦山所能獲取到的數(shù)據(jù)情況，本文選取兩大影響層次，分別為地形地貌層次、地質(zhì)層次。地形地貌層次有坡度、坡向、坡形、高程；地質(zhì)層次有巖性、斷層。

由于本試驗(yàn)采取加權(quán)信息量法，所以在選取因子時(shí)需利用層次分析法確定權(quán)重，并賦予到最后的信息量計(jì)算中。

2.2 層次分析法權(quán)重的確定

層次分析法是19世紀(jì)由Seaty提出的算法，該算法將多目標(biāo)決策問(wèn)題作為一個(gè)系統(tǒng)，將目標(biāo)分為多個(gè)目的和準(zhǔn)則，通過(guò)模糊的定性定量方法將其權(quán)重進(jìn)行計(jì)算［15］。本算法提供3個(gè)層次，分別是目標(biāo)層、決策層和方案層，目標(biāo)層在本次滑坡災(zāi)害風(fēng)險(xiǎn)評(píng)級(jí)里是進(jìn)行滑坡風(fēng)險(xiǎn)的一個(gè)預(yù)測(cè)目標(biāo)，而方案層是所選因子，決策層將所有因子分為地形地貌和地質(zhì)構(gòu)造2個(gè)層次。層次分析法的第二步是將所有不定性不定量數(shù)據(jù)進(jìn)行重要程度比較，也稱專家打分。Seaty將打分表分為幾個(gè)重要程度，如表1所示。

地形地貌層次，參照其他學(xué)者在滑坡案例分析方面選取權(quán)重大小并予以賦值［16］；地質(zhì)層次巖性和斷層的重要程度依舊參考學(xué)者在做地質(zhì)地貌層次打分表進(jìn)行打分［17］。層次分析法在專家打分結(jié)束后，兩兩因子按表1重要程度進(jìn)行比較后便形成層次分析法判斷矩陣，如表2～表4所示。

計(jì)算得出的權(quán)重由于采用專家評(píng)審這一主觀判斷，且隨著層次分析法的判斷矩陣階數(shù)的增大，計(jì)算得出的最大特征向量比矩陣階數(shù)會(huì)大出許多，此時(shí)矩陣的不一致性會(huì)愈發(fā)的加重。為此本文對(duì)采取的打分權(quán)重進(jìn)行一致性分析［18］，從而達(dá)到數(shù)據(jù)的客觀度。層次分析法在進(jìn)行計(jì)算矩陣的特征向量的同時(shí)，也算出各個(gè)層次特征向量所對(duì)應(yīng)的最大特征值λmax，而利用最大特征向量與矩陣階數(shù)的差值（一致性指標(biāo)CI）大小可以衡量矩陣的不一致性。

式中，CI為一致性指標(biāo)；n為層次分析法判斷矩陣階數(shù)。

為防止矩陣因隨機(jī)性而影響一致性指標(biāo)，故結(jié)合隨機(jī)一致性指標(biāo)RI表進(jìn)行參考，一致性指標(biāo)隨矩陣階數(shù)的不確定增大而數(shù)值發(fā)生相應(yīng)的增加（表5）。

結(jié)合一致性指標(biāo)和隨機(jī)引發(fā)矩陣發(fā)生偏差的指標(biāo)，二者之間進(jìn)行比值，便得出我們最終的檢驗(yàn)系數(shù)CR。若CR＜0.1則認(rèn)為該權(quán)重通過(guò)一致性分析，是符合客觀評(píng)價(jià)標(biāo)準(zhǔn)的，該數(shù)據(jù)權(quán)重可被使用；若CR＞0.1則認(rèn)為不通過(guò)一致分析，則需重新取值。

水安全是水資源、水環(huán)境、水生態(tài)、水工程和供水安全五個(gè)方面的綜合效應(yīng)，這五者之間并不是相互獨(dú)立的，如水資源的多寡，與水生態(tài)、供水之間有很強(qiáng)的關(guān)聯(lián)性，水工程、水環(huán)境與供水安全有很強(qiáng)的關(guān)聯(lián)性。

本文2個(gè)方案層，一個(gè)決策層進(jìn)行權(quán)重賦值，其一致性比率分別為：地形地貌因子CR1為0、地質(zhì)構(gòu)造因子CR2為0，總決策層次CR3為0，故權(quán)重通過(guò)一致性檢驗(yàn)，具有可信度。

將得到的各因子權(quán)重按照決策層所占重要程度進(jìn)行歸一化處理如表6所示。

3 加權(quán)信息量法預(yù)測(cè)滑坡

3.1 信息量算法原理及滑坡像元統(tǒng)計(jì)

信息量法是基于信息論方法而建立起來(lái)的一種統(tǒng)計(jì)學(xué)方法，它不僅可以反映各種成災(zāi)要素的相對(duì)敏感程度，也可反映特定成災(zāi)要素中不同要素區(qū)間的致災(zāi)貢獻(xiàn)大?。ㄒ猿蔀?zāi)要素為評(píng)價(jià)因子）［19］。其算法原理：

式中，W為研究區(qū)域內(nèi)單位面積；W0為已經(jīng)發(fā)生滑坡地質(zhì)災(zāi)害的單元面積之和；Sn為具有相同因子x1，x2，x3，…，xn組合的單元總面積；S0為具有相同因子x1，x2，x3，…，xn組合的單元中發(fā)生滑坡地質(zhì)災(zāi)害的總面積。

本文基于大孤山露天礦坑的數(shù)字高程模型（DEM）進(jìn)行像元?jiǎng)澐?，選取DEM柵格精度為0.1 m。為了在地圖中顯示完整，DEM數(shù)據(jù)最大高度與實(shí)際礦坑高程不符，該DEM數(shù)據(jù)是基于實(shí)際高程的500 m以上進(jìn)行建立的。由于所有因子均由大孤山露天礦坑DEM數(shù)據(jù)提取，故本文以DEM數(shù)據(jù)為柵格像元?jiǎng)澐忠罁?jù)，對(duì)DEM柵格進(jìn)行像元統(tǒng)計(jì)。在SuperMap iDesktop中具有柵格統(tǒng)計(jì)功能，統(tǒng)計(jì)得出大孤山露天礦坑總柵格單元為1 676 499個(gè)，根據(jù)現(xiàn)場(chǎng)提供的滑坡位置，進(jìn)行滑坡面的繪制（圖3）。

選取圖3西北幫3個(gè)區(qū)域?yàn)榛挛恢眠M(jìn)行統(tǒng)計(jì)得出，經(jīng)統(tǒng)計(jì)得出滑坡區(qū)域柵格值為12 292個(gè)。

3.2 加權(quán)信息量因子數(shù)據(jù)的提取

基于SuperMap iDesktop自帶的GIS空間分析功能，將大孤山露天礦坑的DEM數(shù)據(jù)進(jìn)行坡度、坡向、坡形、高程的提取，并對(duì)柵格進(jìn)行單元值的分級(jí)。

3.2.1 地形地貌因子

坡向：坡向不同，巖體受陽(yáng)光光照強(qiáng)度不同，土質(zhì)的松散程度也受到相應(yīng)的變化，坡向同樣影響著邊坡失穩(wěn)。坡向提取與坡度提取方法同理，均采用GIS自帶的空間分析功能。以22.5°為間距，進(jìn)行坡向等級(jí)劃分，共劃分16個(gè)等級(jí)。經(jīng)統(tǒng)計(jì)得出90°～135°的坡向發(fā)生滑坡面積最多，占滑坡總面積34%。

坡形：通過(guò)對(duì)DEM表面曲率來(lái)獲取坡形，選取曲率大于0的值為直線型或凸形坡，小于0為凹形坡；經(jīng)統(tǒng)計(jì)得出凸形坡與凹形坡所占滑坡面積數(shù)量相同，二者均勻分布。

高程：選取每100 m為1個(gè)等級(jí)進(jìn)行劃分；共劃分6個(gè)等級(jí)，經(jīng)統(tǒng)計(jì)得出-300～-200 m（實(shí)際高程）發(fā)生滑坡災(zāi)害面積最多，占比為57.13%。

3.2.2 地質(zhì)因子

斷層：距離斷層的大小同樣影響著邊坡的失穩(wěn)程度，由于大孤山露天礦存在斷層數(shù)量較多，故將斷層因子考慮在內(nèi)。以斷層矢量線段為中心，在該中心擴(kuò)展20 m和50 m的圓形區(qū)域（緩沖區(qū)）。斷層分別構(gòu)建了20 m和50 m的緩沖區(qū)，以緩沖區(qū)域進(jìn)行統(tǒng)計(jì)得出距斷層20 m的緩沖區(qū)，發(fā)生滑坡災(zāi)害的面積占比為16.31%；距斷層50 m發(fā)生滑坡災(zāi)害面積占比為26.44%。

巖性：通過(guò)大孤山露天礦地質(zhì)報(bào)告來(lái)進(jìn)行巖性劃分，以不同巖性賦予不同柵格值進(jìn)行云圖的繪制。本文巖性劃分為千枚巖、片麻狀糜棱巖、花崗斑巖、綠泥石英片巖、礦體、閃長(zhǎng)玢巖。其中滑坡災(zāi)害區(qū)域西北幫片麻狀糜棱巖占比最高，為70.33%。

由此便得到了坡度、坡向、坡形、高程、斷層、巖性各自對(duì)應(yīng)的柵格圖。經(jīng)信息量統(tǒng)計(jì)結(jié)果得出，信息量最大值（絕對(duì)值）為1.565 5，所對(duì)應(yīng)因子為距斷層50 m；最小值（絕對(duì)值）為凹形坡，其信息量值為0.033 4。

3.3 加權(quán)信息量法的集成

利用的信息量計(jì)算公式（3）進(jìn)行信息量統(tǒng)計(jì)，再結(jié)合上述層次分析法所確定的各因子的權(quán)重，待獲得上述數(shù)據(jù)后，根據(jù)式（4）可得知最終的加權(quán)信息量［20］，并統(tǒng)計(jì)如表7所示。

注：A1為地形地貌因子，A2為地質(zhì)因子，信息量為負(fù)值時(shí)表示該因子與滑坡發(fā)生與否成負(fù)相關(guān)。

式中，Wi為該因子下層次分析法所算得的權(quán)重；I(y,xi)為該因子下的信息量的值。

4 大孤山露天礦滑坡風(fēng)險(xiǎn)云圖

根據(jù)計(jì)算得出的加權(quán)信息量的值，按照對(duì)應(yīng)柵格所相應(yīng)的加權(quán)信息量賦值，賦值采用GIS軟件的柵格重分級(jí)對(duì)相應(yīng)區(qū)間進(jìn)行相應(yīng)賦值。由此就得到了相應(yīng)因子的加權(quán)信息量的云圖。圖4對(duì)坡度進(jìn)行信息量加權(quán)的柵格進(jìn)行分級(jí)并予以賦值。

在進(jìn)行重分級(jí)后的柵格，便是所求得加權(quán)信息量圖。同理可得坡向、坡形以及其他因素的加權(quán)信息量圖，如圖5所示。

將得到的所有信息量云圖進(jìn)行加權(quán)計(jì)算，由于已經(jīng)賦予了相應(yīng)的層次分析法的權(quán)重，所以只需按照簡(jiǎn)單的線性相加即可。這樣便得到最終的滑坡風(fēng)險(xiǎn)圖。經(jīng)柵格統(tǒng)計(jì)得知信息量的取值范圍為-0.432 96～0.579 44之間，采用自然間斷法分為5級(jí)：低風(fēng)險(xiǎn)級(jí)、次低險(xiǎn)級(jí)、中風(fēng)險(xiǎn)級(jí)、次高、高風(fēng)險(xiǎn)，得出大孤山露天礦山最終滑坡風(fēng)險(xiǎn)圖，如圖6所示。

5 結(jié)論

（1）本次加權(quán)信息量法采取2個(gè)決策層：地質(zhì)層、地形層進(jìn)行分析，通過(guò)GIS空間分析和其他工具手段提取影響災(zāi)害因子，并將其運(yùn)用在礦山，對(duì)礦山災(zāi)害風(fēng)險(xiǎn)評(píng)價(jià)和監(jiān)測(cè)預(yù)警都有著指導(dǎo)意義。

（2）經(jīng)信息量統(tǒng)計(jì)，影響礦坑發(fā)生滑坡最高的因子是坡向在90°～115°因子和距斷層50 m以內(nèi)的因子。

（3）通過(guò)對(duì)大孤山露天礦坑滑坡風(fēng)險(xiǎn)概率云圖分析高風(fēng)險(xiǎn)區(qū)主要在大孤山露天礦西北幫分布，這一點(diǎn)與現(xiàn)場(chǎng)實(shí)際考察較為相符。礦山高風(fēng)險(xiǎn)區(qū)占總礦區(qū)面積比例為7.2%。