李 璐,徐根祺,楊 倩,王艷娥,趙正健
(1.西安思源學院 理工學院,西安 710038;2.西安交通工程學院 機械與電氣工程學院,西安 710030)
滑坡災害的發(fā)生是在自然演變或人為因素的影響下,一種復雜的非線性動力學演化過程[1]?;聻暮τ捎谒旧砀哳l發(fā)生、分布區(qū)域廣泛及破壞力極強,對山區(qū)人民生命、財產有極大的威脅,對防災減災工作提出嚴峻的考驗[2]。近年來國家相關部門和單位陸續(xù)出臺相關政策,滑坡等地質災害的課題也成為熱門課題,相關學者對其研究也取得不錯的成效。滑坡早期預報可以有效減少災害的損失,影響滑坡發(fā)生的主要條件[3-4]:巖土條件、地質構造及地面地形的條件。近年來對滑坡的研究主要包括:1)通過地面調查結合遙感技術觀察并分析滑坡全域地形地貌,分析成災機理[5];2)對滑坡易發(fā)區(qū)域地質構造分析滑坡形成的特征及成因,對滑坡滑動面失穩(wěn)、變形及位移建模[6-7];3)地面外界因素如暴雨、坡體、地震等綜合影響,進行滑坡災害的分析及建模[8-9]。
隨著機器學習理論的發(fā)展,非線性模型被廣泛應用在滑坡災害預測的理論中[10],趙曉萌[11]等從降雨量特征結合機器學習相關知識進行預測模型的建立;趙彬如[12]等從水文、氣象閾值進行滑坡預測的研究,預測效果較好,但只適用于降雨型的滑坡災害預測;胡欣等[13]將SVM-BP模型應用于降雨型的滑坡災害安全評價的模型研究中,隨后馬娟、楊宗佶等[14-15]將多參數預警模式應用于滑坡預測中,同時李麗敏等[16-17]也將多影響因子作為滑坡位移預測模型的輸入,針對滑坡災害建模預報災害的發(fā)生。針對滑坡體災害概率預測及預警的相關研究目前存在以下不足[18]:地域成災機理復雜,成災因子單一,導致預測預警準確率低;預測模型容易陷入局部最優(yōu)及模型本身的穩(wěn)定性及適應性不足。本文針對這些問題分析滑坡全域地形地貌成災機理并篩選因子,構造滑坡災害預測模型,引入核極限學習機(KELM,kernel extreme learning machine)應用于滑坡預測預報中,改善極限學習機(ELM,extreme learning machine)奇點容易產生及隱含層數量確定問題,增加模型本身的穩(wěn)定性及適用性,并使用高斯核函數(RBF,radical basis function)作為核函數進行模型建立。以研究區(qū)山陽縣山區(qū)6種影響因子作為滑坡預測模型的輸入,通過極限梯度提升樹算法Xgboost(Xtreme gradient boosting)確定KELM模型中的參數懲罰系數及寬度參數,解決梯度決策樹算法容易局部最優(yōu)的問題,提高滑坡預測的精度。最后通過與其他預測算法進行比對,AUC均值、Precision、Accuracy及Recall都明顯提升,體現出超參數優(yōu)化后的模型在滑坡預測中較高的預測能力,給地質災害研究帶來活力及新思路。
ELM是一種單隱層的前饋神經網絡,對于線性問題計算速度快且擬合能力也較好,但由于其隱藏節(jié)點存在隨機映射特點,所以訓練的模型穩(wěn)定性和泛化性能力較差[19-20]。KELM在其基礎上將核函數引入解決隨機映射問題。針對于組不同的樣本,ELM隱含層節(jié)點數用表示,激勵函數用g(·)表示:
Hβ=U
(1)
(2)
H為隱含層的輸出矩陣,β為輸出權值矩陣,U為目標輸出矩陣。
ELM學習過程中為提升穩(wěn)定性及泛化能力,引入正則化系數C利用最小二乘法求解最優(yōu)值[21]。
(3)
H+為H的廣義逆。針對ELM解決多分類問題,依據KKT優(yōu)化拉格朗日函數得到求解方案。當特征映射h(x)未知的情況得出ELM的輸出函數:
(4)
對角矩陣用I表示,正則化系數用C表示,輸出向量用U表示。依據Mercer’s條件,用核矩陣ΩELM代替隨機矩陣HHT得出:
K(xi,xj)=h(xi)h(xj)=ΩELM(i,j)
(5)
得到KELM的輸出模型為:
(6)
KELM中核函數會直接影響模型本身的性能。圖1為RBF高斯核函數曲線圖,x為分布中心偏移程度,sigma為分布的寬窄程度,當樣本值接近分布中心x,對核函數值影響較大,反而影響較小,這與內積原理相似,可以通過距離衡量樣本的相似性,高斯徑向基核函數因其是局部核函數,所以局部學習能力強且只受相離較近樣本點的影響,所以其對于一定范圍內距離的樣本可以在特征空間種線性分離預測較準確。文中選取RBF為核函數:
(7)
圖1 RBF高斯核函數
核極限學習機在ELM的基礎上引入了核函數,通過非線性映射到高維特征空間的方式將線性不可分的問題進行劃分,進而提高了ELM性能。但由于核函數的引入,使得KELM算法對參數的選擇非常敏感,所以引入Xgboost算法對KELM正則化系數C和核函數進行尋優(yōu),減少人為調參的復雜度,提高模型的準確性。
Xgboost[22-23]優(yōu)化了梯度提升決策樹,將原來針對錯誤樣本分類中權值分配的不斷迭代,改變?yōu)樨澬牟呗?,訓練最佳方向是損失函數梯度下降的方向,最后通過加權求和的方法得出。并使用增加正則化優(yōu)化了損失函數求解的方法,模型復雜度得到一定降低并防止了過擬合現象,模型精度得到一定程度的提升。
損失函數增加正則化,求解最優(yōu)模型使用結構風險最小的思想[24]設置目標函數來尋優(yōu)迭代模型,如公式(8):
(8)
i為樣本索引符號,I為樣本總量,l為決策樹的葉子節(jié)點數目,Yi為訓練樣本實際值。根據殘差公式fn(Xi)=fn-1(Xi)+hn(Xi),可得出fN(Xi)決策樹訓練(N=1,2,…,N)積累的殘差及Ω(hn(X))針對Ω(h(X))在n次訓練葉子節(jié)點得分,葉子節(jié)點得分ω用L2正則化的表示:
(9)
使用泰勒展開公式進行求解得到二次項,如式(10):
(10)
并規(guī)定:
pi= ?fn-1 (Xi )l(yi,fn-1(Xi))
(11)
最終得出目標函數如(12)所示:
(12)
用該目標函數進行模型最優(yōu)解的求解,可以自己定義損失函數,對于模型構建的靈活性有極大的提升。由于Xgboost基于思想通過不斷最小化目標函數迭代生成決策樹,其預測偏差得到不斷降低,綜合兩方面因素,Xgboost泛化誤差得到整體降低,模型精度勢必得到一定程度的提升。
本文選取裂縫位移ΔX、岸坡水文地質條件H、土壤含水率D、土壓力ΔF、斜坡傾角θ、及降雨量R6個滑坡誘發(fā)條件作為滑坡預測影響因子,監(jiān)測數據通過數據預處理作為模型輸入數據,通過Xgboost優(yōu)化超參數C和并訓練出 Xgboost-KELM模型,分別與GA及GC優(yōu)化的KELM模型比對,最后通過驗證集驗證模型預測的結果。文中具體預測路線如圖2所示。
樣本數據使用陜西省山陽縣2018年3月到2020年3月的10個監(jiān)測點數據作為數據集,樣本集分為80%測試集和20%的兩個驗證集。模型輸入數據:裂縫位移ΔX、岸坡水文地質條件H、土壤含水率D、土壓力ΔF、斜坡傾角θ、及降雨量R6個滑坡體影響因子。預報模型訓練及驗證數集總共1 280組數據集。
表1 監(jiān)測數據
監(jiān)測數據來源于不同的采集傳感器,由于環(huán)境的影響會出現一些缺失、離群或維度不統(tǒng)一的數據,對于模型的建立有極大的消極影響,因此需要對監(jiān)測數據進行預處理。
2.3.1 異常值處理
監(jiān)測數據中存在一部分偏離傳感器本身范圍的值或者偏離觀測值較大的值,如果不處理會影響數據本身預測的準確性,如果距離達到5倍或者相距均值距離≥3倍標準差的數據為離群點。
2.3.2 缺失值的處理
監(jiān)測數據通過多傳感器傳輸,傳輸過程中經常會出現遺漏或者離群點情況,會損失有效信息,導致屬性值確實。按照屬性因素方法進行統(tǒng)計得出缺失率q,本文劃分兩種類別數據的缺失值,如表2所示。
表2 數據缺失值
2.3.3 數據歸一化
采集的監(jiān)測數據種類及數量都較大,多傳感器數據量綱不同有較大的差異,原始數據直接建模對于預測的準確性影響極大,歸一化處理公式如下:
(13)
式中,R為某因素歸一化處理后的數據Rmin,和Rmax為某因素數據中最小值及最大值。圖3為選取4個傳感器中部分監(jiān)測數據數據歸一化前后的數據分布圖,圖中可以看出歸一化前的數據跨度分布比較大,歸一化后的數據在[0~1]量綱內,避免了數據本身量綱問題對預測模型的影響。
圖3 監(jiān)測數據歸一化分布
選取高斯徑向基RBF 核函數,根據(6)和(7)兩式可知需要優(yōu)化的參數有正則化系數C和高斯徑向基RBF 核函數中的核參數σ,其中C∈(0,a],σ∈(0,b]。正則化系數C對模型方差、偏差、訓練誤差及測試誤差都有較大的影響,具體影響如表3所示。核參數σ主要影響模型的擬合程度,當σ>0且較小時容易出現過擬合,相反σ>0較大時容易出現欠擬合。Xgboost優(yōu)化的核極限學習機就是利用 Xgboost優(yōu)化算法對 KELM 中的參數進行擇優(yōu)選取,從而提升 KELM 的性能,提高分類準確度,限制模型復雜度,緩解模型過擬合問題。為了避免多次枚舉造成運算量過大,采用貪心算法尋求最優(yōu)數結構,當Gain信息增益達到樹的深度限制或Gain<0數停止分割,防止過擬合的前提下達到速度快擬合效果好。
表3 正則化參數C對模型指標影響
尋優(yōu)具體流程如下:
STEP 1:數據標準化處理,劃分訓練樣本;
STEP 2:初始化Xgboost算法參數;
STEP 3:Xgboost損失函數使用KELM模型中均方誤差(MSE,mean-square error)代替,選擇Xgboost目標函數下降最大點作為最佳切分點,對RBF 核函數參數σ及正則化系數C進行初始化;
STEP 4:將樣本特征順序排列,列出所有劃分特征、特征值及score分值,根據TOP1分裂子樹,同時計算分割的葉子節(jié)點的權重向量及信息增益;
STEP 5:判斷是否達到數的深度限值或增益小于0,更新并保存最終葉子節(jié)點的權重值與增益值;
STEP 6:判斷最大迭代次數的條件是否達到,如果滿足條件則確定此時的σ及C參數為最優(yōu)參數,基于最優(yōu)參數構建Xgboost-KELM最優(yōu)模型;不滿足條件則返回執(zhí)行STEP 4。
Xgboost參數設置:
1)通用參數
booster基學習器:gbtree(樹模型);
多線程:nthread;
迭代次數nround:1 000。
2)Booster 初始參數
eta:0.1,通過調整學習率調整模型最佳收斂速度;
min_child_weight:0,設置最小葉子節(jié)點樣本的權重和避免出現過擬合現象;
max_depth:8,通過改變樹的最大深度控制子節(jié)點分裂,避免出現過擬合現象;
gamma:節(jié)點分裂的依據,后損失函數下降值大于該值分裂;
subsample:0.7,對樹隨機行比例劃分采樣控制;
n_estimator:120,控制最大樹的數量,防止數量大過擬合,數量小欠擬合;
lambda:1,L2正則化項;
Alpha:1,L1正則化項。
3)目標參數
objective:multi:softmax多分類問題;
eval_metric:[error,auc,mse],回歸模型。
使用訓練集進行訓練Xgboost模型,當Xgboost找尋最優(yōu)的分裂節(jié)點時,可以基于KELM損失函數迭代確定Xgboost最佳參數。圖4為模型篩選最佳參數過程曲線圖,a、b表示最大樹數量與分類準確率關系,a為整個范圍最大樹數量分類過程,最大樹數量n_estimator范圍[1~1 000],變化曲線先逐漸增加,增加至92棵樹時上升緩慢,之后提升不明顯漸漸趨于穩(wěn)定;b為90~150棵樹分類過程,隨之n_estimator數量增加,會存在不到1%的下降趨勢,當數量達到102棵樹時分類準確率達到最高86.34%,取該數量為模型n_estimator最佳參數。c為樹最大深度分類變化過程,max_depth范圍為[1~10],迭代過程可以看出深度為6時分類準確度最高達86.42%,max_depth最佳參數值取6;d為最小葉子權重變化過程,葉子節(jié)點的權重小于min_child_weight則停止拆分樹,min_child_weight范圍為[0~10],曲線變化過程可以得出權重為2時分類準確率達到最高88.25%,min_child_weight最佳值取2。根據訓練集訓練結果得出Xgboost最終參數如表4所示。
圖4 不同參數分類準確率
表4 Xgboost模型最終參數
為了獲取KELM最優(yōu)核參數σ及正則化C參數最佳組合,將測試集隨機分為4組作為訓練樣本建立模型,參數的選值先從一個比較大的區(qū)間范圍進行搜索,4組測試集驗證后得出小范圍[500,1 000]×[0.2,0.3],之后進行多次迭代訓練,當均方誤差達到最小值時,取此時參數核參數σ及正則化C為最佳組合參數。圖5為4組測試集訓練過程的均方誤差迭代曲線,數據集4收斂速度最慢,數據集2、3均方誤差較低但收斂速度在迭代次數12出現飽和,而數據集3均方誤差及收斂速度相比其他數據集表現最佳,均方誤差最低達到并穩(wěn)定于1.187×10-3。且在7次迭代收斂飽和。所以選取數據集3訓練參數為模型最佳參數,σ=0.262 4,C=703.24。
圖5 不同測試集的均方誤差
為驗證模型的穩(wěn)定性及適應能力,通過模型在新樣本集中的適應度、準確性及方差偏差驗證。方差表示模型每次預期結果與實際結果的誤差的穩(wěn)定情況;偏差值表示每次預期結果與實際結果的偏差。模型誤差包括方差、偏差及其他無法避免誤差,圖6為偏差及方差示意圖。預測模型最佳選擇順序:1)方差小,偏差??;2)方差小,偏差大;3)方差大,偏差?。?)方差大,偏差大。
圖6 偏差及方差示意圖
精確率:在預期的正樣本中實際結果也為正樣本的占比。
precision=TP/(TP+FP)
(14)
準確率:準確率表示所有的預測樣本中,預測正確的占比。
Accuracy=(TP+TN)/(TP+FP+FN+TN)
(15)
召回率:預測結果準確的正樣本占所有正樣本的比例。
recall=TP/(TP+FN)
(16)
AUC:通過計算ROC曲線與坐標軸圍成的面積得到,介于[0.5,~1]之間,預測的真實性取決與AUC值接近1的程度,靠近1真實性高反之則反。
(17)
真正率:預期正樣本數/實際正樣本數。
(18)
假正率:預期為正的負樣本數/實際負樣本數。
(19)
實驗采用KELM作為滑坡災害預測模型,并用Xgboost優(yōu)化算法中的超參數尋優(yōu)。使用同一個驗證集驗證GC及GA優(yōu)化KELM模型預測效果比對。表5為在驗證集1中128個數據中各模型輸入相同數據得出的預測結果的混淆矩陣對比表。
表5 不同模型混淆矩陣對比表
圖7為3種模型同一評價指標的對比圖,柱狀圖的差異可以對比得出各個預測模型對應評價指標的好壞,從而反映預測模型性能的優(yōu)劣。4個評估參數中,Xgboost優(yōu)化KELM均明顯優(yōu)于GA和GC優(yōu)化的模型。本文的Xgboost-KELM模型AUC均值為0.985,相比GC優(yōu)化高約3個百分點,比GA優(yōu)化高6個百分點。其他指標Precision、Accuracy及Recall高于另外兩個模型百分比[1.7~7]范圍之間。實驗結果說明Xgboost-KELM模型具有較好的預測效果,在滑坡災害預測中有較好的預測能力。
圖7 評價指標對比圖
圖8為驗證集2中打亂隨機抽取的128個監(jiān)測數據使用Xgboost優(yōu)化KELM模型后的實際發(fā)生概率與預測發(fā)生概率比對圖。圖中實際值與預測值基本吻合,擬合情況較好。極限的幾個數據27、38、89及111發(fā)生概率存在一些差異,但是對應風險等級都屬于同等級風險。準確率達到98%。引入Xgboost對KELM參數值進行優(yōu)化,實際預測精度較為理想。并使用驗證集2進行各模型穩(wěn)定性計算,Xgboost-KELM有較小的方差且偏差也較小,屬于最穩(wěn)定的“方差小、偏差小”模型,穩(wěn)定性能最強,而GC-KELM介于“方差大、偏差小”與“方差大、偏差大”之間,模型不穩(wěn)定,GA-KELM屬于“方差大、偏差大”,模型最不穩(wěn)定。
圖8 Xgboost優(yōu)化模型預測結果
本文針對山陽縣研究區(qū)域的數據使用基于Xgboost優(yōu)化KELM模型建立滑坡災害預測模型,通過仿真研究分析滑坡影響因子與發(fā)生概率之間關系,并與GA-KELM及GC-KELM建模結果進行比較。
1)選用RBF高斯核函數的極限學習機模型較好地解決了ELM適用性及穩(wěn)定性不佳的問題;
2)使用KELM模型中均方誤差MSE作為損失函數,訓練模型并選擇目標函數下降最大點作為最佳切分點,確定最佳參數;并使用Xgboost尋優(yōu)算法對核函數中的正則化系數C和核函數σ尋優(yōu),通過4組測試集的均方誤差迭代曲線得出最佳超參數,建立Xgboost-KELM模型,并與GA及GC優(yōu)化KELM建模進行比較;
3)通過幾種模型的樣本集1對比驗證,結果表明Xgboost-KELM具有較高的Precision、Accuracy及Recall和AUC值,同時使用新樣本集2驗證穩(wěn)定性及泛化能力,結果表明該模型穩(wěn)定性較好,進一步證明該模型的應用能有效提高滑坡災害的預報概率,對滑坡災害提前預警,降低自然災害造成的損失具有重要意義。