林逸暉,李廣濤,楊天雨,喬登攀,王 俊,張 希,趙懷軍

(1.昆明理工大學(xué) 國土資源工程學(xué)院,云南 昆明 650093;2.涼山礦業(yè)股份有限公司,四川 涼山 615141)

0 引言

近年來,隨著現(xiàn)代工業(yè)對礦產(chǎn)資源的需求量日益增加,我國露天礦山的開采規(guī)模不斷擴(kuò)大,開采深度越來越大,高陡邊坡數(shù)量越來越多,邊坡災(zāi)害發(fā)生頻次逐年上升。邊坡災(zāi)害會對人民生命財(cái)產(chǎn)安全造成嚴(yán)重威脅,并使生態(tài)環(huán)境遭受巨大破壞[1-2]。因此,為減少邊坡災(zāi)害的發(fā)生,開展邊坡穩(wěn)定性預(yù)測研究具有重要的現(xiàn)實(shí)意義。

露天礦山邊坡穩(wěn)定性分析方法主要有定性分析法、定量分析法和數(shù)值分析法等[3]。隨著計(jì)算機(jī)技術(shù)的發(fā)展,機(jī)器學(xué)習(xí)算法被引入邊坡穩(wěn)定性預(yù)測研究中。目前用于邊坡穩(wěn)定性預(yù)測的機(jī)器學(xué)習(xí)算法有粗糙集理論、模糊邏輯理論、支持向量機(jī)、蟻群算法、貝葉斯算法、信息量法、多元邏輯回歸、決策樹法等[4]。瞿婧晶等[5]采用粗糙集屬性約簡和相互關(guān)系矩陣法進(jìn)行綜合分析,建立了邊坡穩(wěn)定性評價體系,并將其應(yīng)用于鎮(zhèn)江地區(qū)下蜀土邊坡的穩(wěn)定性預(yù)測中。張?jiān)蒲鉡6]采用網(wǎng)格搜索算法對支持向量機(jī)模型進(jìn)行參數(shù)優(yōu)化,并根據(jù)邊坡特征參數(shù),采用優(yōu)化后的支持向量機(jī)模型對邊坡穩(wěn)定性進(jìn)行了預(yù)測。趙允坤等[7]利用改進(jìn)粒子群算法(IPSO)對極限學(xué)習(xí)機(jī)的輸入權(quán)值和隱層偏置進(jìn)行了優(yōu)化,建立了基于IPSO-ELM的邊坡穩(wěn)定性預(yù)測模型,實(shí)現(xiàn)了對邊坡穩(wěn)定性的有效預(yù)測。

但機(jī)器學(xué)習(xí)算法也存在一定缺陷,如:粗糙集理論需結(jié)合其他算法使用,否則難以有效反映不確定性問題;支持向量機(jī)模型對于缺失數(shù)據(jù)以及參數(shù)和核函數(shù)的選擇較敏感,且核函數(shù)及其參數(shù)是根據(jù)經(jīng)驗(yàn)選取的,帶有一定的主觀性;極限學(xué)習(xí)機(jī)存在易過擬合以及可控性差等缺點(diǎn)。

隨機(jī)森林(Random Forest,RF)是由美國科學(xué)家Leo Breiman[8]于2001年提出的,這是一種將Bagging集成學(xué)習(xí)理論與隨機(jī)子空間方法相結(jié)合的機(jī)器學(xué)習(xí)算法。隨機(jī)森林以決策樹為基本分類器,相對于決策樹法,隨機(jī)森林克服了其易出現(xiàn)過擬合的缺點(diǎn),降低了異常值可能帶來的影響,但計(jì)算量較大[9-10]。此外,隨機(jī)森林對高維數(shù)據(jù)分類問題有著良好的可擴(kuò)展性和并行性,并且作為一種由數(shù)據(jù)驅(qū)動的非參數(shù)分類方法,使用隨機(jī)森林算法時無需調(diào)參以及分類先驗(yàn)知識[11]。針對隨機(jī)森林算法計(jì)算性能開銷較大的缺點(diǎn),本文利用主成分分析法(Principal Component Analysis,PCA)進(jìn)行數(shù)據(jù)處理,過濾冗雜數(shù)據(jù),降低隨機(jī)森林算法的計(jì)算量,構(gòu)建基于PCA-RF的邊坡穩(wěn)定性預(yù)測模型,并將其應(yīng)用于露天邊坡的穩(wěn)定性預(yù)測。

1 理論基礎(chǔ)

1.1 主成分分析

主成分分析是利用降維思想,在信息損失較少的前提下,將原來具有一定相關(guān)性的指標(biāo)重新組合成一組新的信息互不重疊的綜合指標(biāo),來反映大部分原指標(biāo)中所攜帶的信息量的多元統(tǒng)計(jì)方法[12]。主成分分析計(jì)算步驟敘述如下。

(1)建立m×n矩陣,m為樣本數(shù)量,n為每個樣本的指標(biāo)個數(shù)。

(1)

(2)將原始數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化,生成標(biāo)準(zhǔn)化矩陣。

(2)

(3)進(jìn)行KMO(Kaiser-Meyer-Olkin)和巴特利特球形檢驗(yàn)。若KMO值大于0.5,表明變量間的相關(guān)程度差異較小,數(shù)據(jù)很適合作因子分析[13];若巴特利特球形檢驗(yàn)的結(jié)果小于0.05,球形假設(shè)被拒絕,表明原始變量之間存在相關(guān)性,數(shù)據(jù)適合作因子分析[14]。

(4)根據(jù)標(biāo)準(zhǔn)化矩陣建立協(xié)方差矩陣。

(3)

式中,i=1,2,…,n,j=1,2,…,n。

協(xié)方差矩陣由標(biāo)準(zhǔn)化矩陣的兩兩列變量的協(xié)方差組成,因此cov(zi,zj)為列變量zi與zj的相關(guān)系數(shù),即協(xié)方差矩陣為相關(guān)系數(shù)矩陣。

(5)計(jì)算協(xié)方差矩陣R的非負(fù)特征根λi(i=1,2,…,p)。

|λE-R|=0。

(4)

(6)由計(jì)算得到的p個非負(fù)特征根計(jì)算前q個主成分的累計(jì)貢獻(xiàn)率。

(5)

(7)根據(jù)主成分的累計(jì)貢獻(xiàn)率選取主成分個數(shù),并生成新變量指標(biāo)。前q個主成分的累計(jì)貢獻(xiàn)率達(dá)到85%時,即可認(rèn)為這些主成分中包含了絕大部分原指標(biāo)所攜帶的信息量[15]。將生成的新變量指標(biāo)表示原變量指標(biāo):

(6)

式中,x1,x2,…,xn為原變量指標(biāo),y1,y2,…,yp為新變量指標(biāo),cij與yi互不相關(guān)且cij滿足ci12+ci22+…+cin2=1。由此完成了數(shù)據(jù)降維,在信息損失較少的前提下減少了變量數(shù)。

1.2 隨機(jī)森林模型

隨機(jī)森林{h(X,θk),k=1,2,…,K}是以K個決策樹模型作為基本分類器,通過集成學(xué)習(xí)形成的組合分類器[16]。隨機(jī)森林算法以待分類樣本作為輸入,基于決策樹模型,采用自上而下的遞歸方式,從樹的根節(jié)點(diǎn)開始,在其內(nèi)部節(jié)點(diǎn)上進(jìn)行屬性值的測試比較,然后按照給定的屬性值確定對應(yīng)分支,在決策樹的葉子節(jié)點(diǎn)得到結(jié)果,最后根據(jù)每個決策樹模型的分類結(jié)果投票決定最終輸出的分類結(jié)果。

隨機(jī)森林中的{θk,k=1,2,…,K}為隨機(jī)變量序列,其體現(xiàn)在兩個方面:

(1)數(shù)據(jù)的隨機(jī)選取。在Bagging集成算法的基礎(chǔ)上,利用Bootstrap方法從原始樣本集中采取有放回抽樣,隨機(jī)創(chuàng)建K個與原樣本集數(shù)據(jù)量相同的子數(shù)據(jù)集{Tk,k=1,2,…,K},其中不同子數(shù)據(jù)集的元素可以重復(fù),同一個子數(shù)據(jù)集中的元素也可以重復(fù)。

(2)待選屬性的隨機(jī)選取。采用隨機(jī)子空間方法,在隨機(jī)森林中的子決策樹的各節(jié)點(diǎn)進(jìn)行分裂時,從所有的待選屬性中等概率隨機(jī)選取一定的屬性[通常選取(log2M+1)個屬性,M為屬性總個數(shù)[17]],再從所選屬性中選取最優(yōu)屬性進(jìn)行節(jié)點(diǎn)分裂。

隨機(jī)森林的分類過程為:

(1)利用Bootstrap方法隨機(jī)選取數(shù)據(jù),生成子數(shù)據(jù)集Tk,將Tk作為訓(xùn)練集構(gòu)建K棵決策樹。

(2)每棵決策樹生長過程中在所有的待選屬性中等概率隨機(jī)選取一定的屬性,并從中選取最優(yōu)屬性對每棵決策樹的各個節(jié)點(diǎn)進(jìn)行分裂,并且讓每棵決策樹均不進(jìn)行剪枝地完整生長。

(3)由完整生長的K棵決策樹構(gòu)成隨機(jī)森林。

(4)輸入測試集樣本進(jìn)行預(yù)測,由每棵決策樹的輸出結(jié)果投票決定隨機(jī)森林輸出的分類結(jié)果(即每棵決策樹輸出結(jié)果的眾數(shù))。

2 PCA-RF預(yù)測模型構(gòu)建及驗(yàn)證

2.1 樣本數(shù)據(jù)

影響邊坡穩(wěn)定性的因素主要有邊坡角度、邊坡高度、巖土體物理力學(xué)性質(zhì)、巖石強(qiáng)度、巖石結(jié)構(gòu)、孔隙水壓、水文條件及地震等[18]。但巖石強(qiáng)度、巖石結(jié)構(gòu)等因素對于邊坡穩(wěn)定性的影響較小,地震等自然災(zāi)害出現(xiàn)的概率較低,樣本較少。因此,為減少隨機(jī)森林算法的計(jì)算量,提高計(jì)算速度,選取巖石容重γ、黏聚力c、內(nèi)摩擦角φ、邊坡角度α、邊坡高度H和孔隙水壓μ等6個主要影響因素進(jìn)行邊坡穩(wěn)定性預(yù)測[19]。本文基于此6項(xiàng)指標(biāo)構(gòu)建隨機(jī)森林邊坡穩(wěn)定性預(yù)測模型,將使用主成分分析降維后的指標(biāo)作為實(shí)際輸入、邊坡穩(wěn)定狀態(tài)“穩(wěn)定”和“失穩(wěn)”作為輸出進(jìn)行邊坡穩(wěn)定性預(yù)測。

本文選取文獻(xiàn)[20-21]中的54組不同且具有明確穩(wěn)定狀態(tài)的邊坡實(shí)例作為研究對象,構(gòu)建PCA-RF邊坡穩(wěn)定性預(yù)測模型。邊坡實(shí)例數(shù)據(jù)見表 1,其中:邊坡狀態(tài)“1”表示邊坡穩(wěn)定,“2”表示邊坡失穩(wěn);穩(wěn)定邊坡為34組,失穩(wěn)邊坡為20組,第55-第58組數(shù)據(jù)為拉拉銅礦東部露天采場終了邊坡數(shù)據(jù),以此作為預(yù)測實(shí)例。

表1 邊坡實(shí)例數(shù)據(jù)

2.2 主成分分析

對表 1中的58組數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理,計(jì)算得到KMO值為0.665,大于0.5;巴特利特球形檢驗(yàn)值為0.000,小于0.05,表明樣本數(shù)據(jù)適合進(jìn)行主成分分析。對表 1中的數(shù)據(jù)進(jìn)行主成分分析,得到各影響因素間的相關(guān)矩陣熱圖(見圖 1)。將相關(guān)性系數(shù)取絕對值后,如在0～0.1,表示無相關(guān);如在0.1～0.3,表示弱相關(guān);如在0.3～0.5,表示中等相關(guān);如在0.5～1.0,表示強(qiáng)相關(guān)[22]。

圖 1 相關(guān)矩陣熱圖

由圖 1可知,某些影響因素間存在較強(qiáng)的相關(guān)性,如巖石容重與邊坡角度、邊坡高度間的相關(guān)性系數(shù)分別為0.535、0.679,均大于0.5,表明具有強(qiáng)相關(guān)性。主成分分析能夠?qū)⑦@些指標(biāo)重新組合成一組新的信息互不重疊的綜合指標(biāo),從而避免相關(guān)性過強(qiáng)可能帶來的多重共線性問題。

為確保選取的指標(biāo)中攜帶原樣本數(shù)據(jù)的絕大部分信息量,需計(jì)算各因素對邊坡穩(wěn)定性影響的貢獻(xiàn)率及累計(jì)貢獻(xiàn)率,結(jié)果見表2。

表2 各影響因素貢獻(xiàn)率

表3 主成分系數(shù)矩陣

由表 2可知,前4項(xiàng)主成分的累計(jì)貢獻(xiàn)率達(dá)到了91.283%,大于85%,表明前4項(xiàng)主成分已攜帶原樣本數(shù)據(jù)的絕大部分信息量,因此選取前4項(xiàng)主成分生成新變量指標(biāo)。各主成分系數(shù)矩陣見表 3。

將各主成分系數(shù)分別除以相應(yīng)的主成分特征根的平方根,得到因子載荷矩陣(見表4)。

表4 因子載荷矩陣

由因子載荷矩陣可得出新變量指標(biāo)的表達(dá)式:

(7)

經(jīng)主成分分析后,在降低變量維數(shù)的同時,還保證了原樣本數(shù)據(jù)攜帶的信息量盡可能少丟失。經(jīng)主成分分析降維后的邊坡數(shù)據(jù)見表5。

表5 經(jīng)主成分分析降維后的邊坡數(shù)據(jù)

表6 PCA-RF模型預(yù)測結(jié)果

2.3 建立邊坡穩(wěn)定性預(yù)測模型

從降維后的前54組邊坡數(shù)據(jù)中隨機(jī)選取45組數(shù)據(jù)作為輸入訓(xùn)練樣本建立PCA-RF邊坡穩(wěn)定性預(yù)測模型,將剩下的9組數(shù)據(jù)作為測試樣本,用來檢驗(yàn)該預(yù)測模型的準(zhǔn)確率。

在使用隨機(jī)森林進(jìn)行預(yù)測之前,需先設(shè)置2個重要參數(shù):隨機(jī)森林中包含的決策樹個數(shù)ntree和每棵樹節(jié)點(diǎn)預(yù)選的變量指標(biāo)個數(shù)mtry。若預(yù)設(shè)的ntree值太小會導(dǎo)致隨機(jī)森林模型訓(xùn)練不夠充分,太大則會導(dǎo)致模型運(yùn)行速度緩慢;若預(yù)設(shè)的mtry值太小可能會導(dǎo)致分類器過擬合,使預(yù)測精度降低,太大則會增加模型的運(yùn)算量。綜合考慮,本文將ntree值設(shè)為500,mtry值設(shè)為4。PCA-RF模型預(yù)測結(jié)果見表 6。

由表 6可知,PCA-RF模型的預(yù)測準(zhǔn)確率達(dá)到了100%。但因隨機(jī)森林生成訓(xùn)練集時采用Bootstrap有放回自助采樣法,輸入訓(xùn)練樣本中有(1-1/N)N的樣本不會被抽取(N為輸入訓(xùn)練樣本中的樣本個數(shù)),當(dāng)N足夠大時,(1-1/N)N將收斂于1/e,約為0.368,即有約36.8%的樣本數(shù)據(jù)不會被抽取[23]。因此,該預(yù)測模型有可能存在誤差,為此需要對該模型進(jìn)行性能評估。

2.4 交叉驗(yàn)證

模型性能的評估方法有留一法和交叉驗(yàn)證法等,本文選用6折交叉驗(yàn)證法對PCA-RF模型進(jìn)行性能評估。6折交叉驗(yàn)證法驗(yàn)證步驟為:

(1)將54組樣本數(shù)據(jù)隨機(jī)分為數(shù)據(jù)量相同的6組數(shù)據(jù)子集D1-D6。

(2)取其中的第i組數(shù)據(jù)子集Di作為測試樣本,將其余5組數(shù)據(jù)子集作為輸入訓(xùn)練樣本,構(gòu)建6個隨機(jī)森林預(yù)測模型。

(3)匯總預(yù)測結(jié)果,開展PCA-RF模型的性能評估。

交叉驗(yàn)證過程中預(yù)測模型的參數(shù)設(shè)置與前文相同,驗(yàn)證結(jié)果見表7。由表7可知,模型的預(yù)測精度較高,達(dá)到了94.44%,表明該預(yù)測模型可以用于邊坡穩(wěn)定性預(yù)測。

表 7 交叉驗(yàn)證結(jié)果

3 工程應(yīng)用

3.1 預(yù)測對象選取

拉拉銅礦礦體主要部分埋藏不深,礦山原設(shè)計(jì)采用露天開采,露天開采區(qū)域分為東部露天采場和西部露天采場,東部露天采場尚未開采至終了境界,其終了邊坡穩(wěn)定狀態(tài)未知,因此本文選取東部露天采場終了邊坡作為預(yù)測對象,現(xiàn)場照片見圖2。該礦東部露天采場終了邊坡東幫設(shè)計(jì)高度為252 m,坡角為45.57°;西幫設(shè)計(jì)高度為216 m,坡角為41.75°;南幫設(shè)計(jì)高度為296 m,坡角為45.89°;北幫設(shè)計(jì)高度為312 m,坡角為30.05°。露天采場區(qū)域巖性主要為片巖,巖石力學(xué)參數(shù)見表 1中第55-第58組數(shù)據(jù)。東部露天采場邊坡現(xiàn)場照片見圖3。

圖2 拉拉銅礦東部露天采場

圖3 拉拉銅礦東部露天采場邊坡

3.2 預(yù)測結(jié)果分析

將經(jīng)主成分分析降維后的拉拉銅礦東部露天采場終了邊坡數(shù)據(jù)作為預(yù)測數(shù)據(jù),輸入已訓(xùn)練好的隨機(jī)森林預(yù)測模型中,輸出的預(yù)測結(jié)果見表8。由表8可知,拉拉銅礦東部露天采場終了邊坡的東幫、西幫、南幫、北幫均處于穩(wěn)定狀態(tài)。

表8 拉拉銅礦東部露天采場終了邊坡穩(wěn)定性預(yù)測結(jié)果

4 結(jié)論

a.本文采用主成分分析對原始邊坡數(shù)據(jù)進(jìn)行降維,減少了隨機(jī)森林算法的計(jì)算量,建立了PCA-RF邊坡穩(wěn)定性預(yù)測模型。

b.邊坡實(shí)例預(yù)測結(jié)果表明,PCA-RF模型預(yù)測準(zhǔn)確率達(dá)到100%,并且該模型經(jīng)6折交叉驗(yàn)證法評估的預(yù)測精度高達(dá)94.44%,表明該模型可以用于邊坡穩(wěn)定性預(yù)測。

c.將PCA-RF模型應(yīng)用于拉拉銅礦東部露天采場邊坡穩(wěn)定性預(yù)測,結(jié)果表明,該礦東部露天采場東幫、西幫、南幫、北幫終了邊坡均處于穩(wěn)定狀態(tài),該預(yù)測結(jié)果可為礦山生產(chǎn)和管理提供參考。