李 璐，張永強，李麗敏，馬 媛，竇婉婷，王 悅

(1.西安思源學院 理工學院，西安 710038；2.銅川職業(yè)技術學院 機電工程學院，陜西 銅川 727031；3.西安工程大學 電子信息學院，西安 710048；4.西安交通大學城市學院 傳播系，西安 710018)

0 引言

泥石流災害的發(fā)生是在自然演變或人為因素的影響下，一種復雜的非線性動力學演化過程。我國山區(qū)較多，泥石流災害是山區(qū)常見的一種自然災害，由于它本身高頻發(fā)生、分布區(qū)域廣泛及破壞力極強，對山區(qū)人民生命、財產(chǎn)有著極大的威脅，對防災減災工作提出嚴峻的考驗。泥石流的早期預報可以有效減少災害的損失，泥石流形成主要有三大條件，分別是地形地貌、松散物源、水源[1]。近年來國家對地質(zhì)災害的防災減災比較重視，陸陸續(xù)續(xù)出臺政策，隨著災害的頻繁發(fā)生，對泥石流災害的研究一直都是熱度較高的課題，相關學者針對泥石流的研究主要有：1)通過災害區(qū)域地面調(diào)查結(jié)合相關遙感技術，觀察并分析泥石流災害全域的地形地貌，從而分析其成災機理[2]；2)對物源動儲量、泥沙補給、流量等影響因素通過力學及流變學的理論建立相關泥石流的運動方程[3-4]；3)通過實時監(jiān)測收集雨量信息，對降雨強度與臨界雨量閾值分析并建立雨量模型[5]；4)通過實時監(jiān)測采集成災因子，對泥石流發(fā)生的概率及等級進行預報，從而達到提前預報預警提示，減少災害重大損失[6]。隨著機器學習理論的不斷發(fā)展，非線性模型也被廣泛應用在泥石流災害預測的理論中，文獻[7]融合泥石流的多個影響因子，通過遺傳規(guī)劃法建立臨界降雨指數(shù)智能預測模型；文獻[8]基于PCA(principal component analysis)篩選泥石流災害成災因子并使用BP(back propagation)神經(jīng)網(wǎng)絡對泥石流發(fā)生的危險性進行預測，此方法選用有效成災因子的方法結(jié)合預測模型極大提升了泥石流危險性的預測，但是使用PCA篩選因子處理非線性關系有一定缺陷。文獻[9]使用混合核函數(shù)改進了KPCA篩選因子算法，預測等級達到一定的提升?；谶@一思想，學者[10-12]將成災因子篩選、模型參數(shù)尋優(yōu)等優(yōu)化模型方式使得預測模型精度更加穩(wěn)健。本文借鑒這一思想分析泥石流全域地形地貌成災機理并篩選因子，構(gòu)造出泥石流災害預測模型。

為進一步提升泥石流預測的精度，本文以陜西省山陽縣中村鎮(zhèn)泥石流全域為研究對象，首先分析災害區(qū)域地形地貌選出成災因子，避免因使用單因子預測導致的精度低問題；其次采用核主成分分析法(KPCA，Kernel principal component analysis)進行高維度影響因子的篩選；另外構(gòu)造最小二乘支持向量機(LSSVM，least squares support vector machines)模型對泥石流發(fā)生概率模型的建立，相對于支持向量機將非線性問題轉(zhuǎn)換為線性問題求解方式極大的簡化，同時使用多算法進行LSSVM中的超參數(shù)優(yōu)化，經(jīng)過優(yōu)化后的模型解決了過早收斂導致陷入局部最優(yōu)的問題以及參數(shù)隨機選取導致的精度不佳問題，通過模型優(yōu)化提高了泥石流預測的精度。最后通過與其他尋優(yōu)預測算法進行比對，比對出最佳預測模型，為泥石流地質(zhì)災害研究帶來活力及新思路。

1 算法理論

1.1 KPCA方法

主成分分析方法[13]是一種無監(jiān)督降維算法，針對線性數(shù)據(jù)效果較好，但是其對于數(shù)據(jù)之間存在非線性關系時PCA降維效果比較差，為了彌補這一缺陷，在PCA計算協(xié)方差矩陣時加入核函數(shù)，用來解決非線性映射問題。KPCA在高維特征空間對原數(shù)據(jù)映射，經(jīng)過PCA對高維數(shù)據(jù)計算特征向量及特征值來確定主成分因子。

(1)

通過特征分解C值得出：

λν=Cν

(2)

當所有特征值λ≥0，v為由φ(x1)，φ(x2)，…，φ(xM)組成的空間，所以式(2)等于：

λ(φ(xk)，vr)=(φ(xk)，Cvr)k=1，2，…，M

(3)

由于vT是φ(x)的線性組合，所以得出：

(4)

將式(1)、(4)代入式(3)，并令Kij=(φ(xi)，φ(xj))i，j=1，2，…，M，代入得出：

Mλrcr=Kcr

(5)

Mλr為特征值，cr為特征向量，當滿足cr>0條件：cp，cp+1，…，cM，進行歸一化后得出：

Mλr(cr，cr)=1

(6)

求得φ(x)在cr特征向量的投影：

r=p，p+1，…，M

(7)

g(x)為φ(x)非線性主元分量，g(x)[g1(x)，g2(x)，…，gI(x)]T為所有投影矢量表示。使用核函數(shù)K(xi，xj)=〈φ(xi)，φ(xj)〉求解g(x)代替空間的點積運算，核函數(shù)變?yōu)椋?/p>

g(x)=(vT，φ(x))=K(xi，xj)

(8)

當φ(x)≠0時，空間樣本變換：

(9)

通過式(8)計算矩陣K，再依據(jù)樣本變換求取特征向量與特性值，最后依據(jù)最大特征值及其對應向量結(jié)合輸入屬性得到主成分。按照式(10)、(11)得出各個成分的貢獻率與累計貢獻率。

(10)

(11)

1.2 LSSVM模型

LSSVM(least squares support vector machines)[14-15]基于SVM將不等式約束轉(zhuǎn)換為等式約束，從而化簡lagrange乘子α求解，對求解QP問題轉(zhuǎn)為進行線性方程組的求解。LSSVM 繼承了 SVM 的泛化能力和魯棒性，但其計算效率優(yōu)于原始的 SVM。給定訓練的數(shù)據(jù)集合(xi，yj)，i=1，2，…，n，分別給出SVM及LSSVM需求解的問題。

SVM不等式約束問題：

s.t.yi[ωT·φ(xi)+b]≥1-ζi，i=1，2，…，n

(12)

LSSVM等式約束問題：

s.t.yi[ωT·φ(xi)+b]=1-ei，i=1，2，…，n

(13)

ξ及e為松弛變量，用于SVM及LSSVM中引入離群點，c及γ為平衡尋找最優(yōu)超平面與偏差量之間最小值，ω為權重向量，b為誤差，φ(·)為映射函數(shù)。

使用Lagrangea方法對式(13)優(yōu)化，轉(zhuǎn)化為單一的參數(shù)，求解α的極限值，構(gòu)造出：

L(ω，b，e，α)=

(14)

其中：αi為拉格朗日乘子。

ω，b，ei，αi分別求導=0：

(15)

φ依據(jù)4個求導的條件可列出線性方程組：

(16)

In為單位矩陣的轉(zhuǎn)置矩陣，E為n維單位矩陣，Φ為核矩陣：

Φij=yiyjφ(xi)Tφ(xj)=

yiyjK(xi，xj)，i，j=1，2，…，n

(17)

解方程(16)，可得出一組α、b，最后得出LSSVM分類表達式為：

(18)

LSSVM的訓練框架如圖1所示，LSSVM算法中的正則化系數(shù)和核函數(shù)參數(shù)需要進行尋優(yōu)防止出現(xiàn)參數(shù)隨機導致的精度不佳問題及過早收斂導致陷入局部最優(yōu)的問題。

圖1 LSSVM模型網(wǎng)絡結(jié)構(gòu)圖

1.3 松鼠覓食算法

松鼠覓食算法(sparrow search algorithm)[16]對于搜索空間中的一些復雜問題搜索能力及精度有明顯優(yōu)勢，松鼠雖不會飛行，但可以通過滑翔的方式來躲避天敵捕食，SSA算法就是模擬其這一行為的過程。松鼠的搜尋過程伴隨其覓食的開始，尋找食物的方式通過其從不同的樹木之間移動來獲取，森林中不同區(qū)域的搜索通過松鼠位置的改變來實現(xiàn)。

假設松鼠的數(shù)量為n，松鼠移動的位置通過矢量來確定，并在邊界范圍內(nèi)隨機初始化其位置。

(19)

FSn，d為第n只老鼠在第d維度上的值，松鼠在森林中的初始位置為：

FSi=FSL+U(0，1)×(FSU-FSL)

(20)

FSU和FSL為松鼠移動的上下界，U(0，1)為隨機數(shù)[0，1]。

食物源的等級通過每一只松鼠位置的適應度表示，計算適應度值并進行升序分類，適應度最小的位置：最佳食物源①山核桃，接下來三只位置正常食物源②橡樹，其他的位置無食物來源③普通樹。

依據(jù)其天敵的出現(xiàn)概率Pdp松鼠更新移動的位置。

滑翔路徑一：②→①

(21)

滑翔路徑二：③→②

(22)

滑翔路徑三：③→①

(23)

季節(jié)的變化會影響松鼠的覓食活動，使用季節(jié)的變換來防止出現(xiàn)算法陷入局部最優(yōu)。

(24)

(25)

(26)

FSi，U和FSi，L為松鼠移動的上下界，Levy為列維分布，有效地全局搜索，來找到距離當前地點最優(yōu)的一個新地點。SSA算法具體步驟如圖2所示。

圖2 SSA算法流程

2 基于KPCA-SSA-LSSVM山區(qū)泥石流災害預測模型

基于KPCA-SSA-LSSVM的山區(qū)泥石流災害發(fā)生預測流程如圖3所示，具體實現(xiàn)步驟如下：

圖3 KPCA-SSA-LSSVM算法流程

1)首先對監(jiān)測的原始數(shù)據(jù)預處理，并使用KPCA核主成分分析法篩選出覆蓋率滿足需求的6個影響因子。

2)數(shù)據(jù)集合理劃分，確定訓練集及測試集百分比。

3)初始化尋優(yōu)參數(shù)及LSSVM參數(shù)。

4)根據(jù)各影響因子建立LSSVM預測模型，并在訓練集訓練最佳適應度函數(shù)。

5)將不同預測模型在測試集上分析對比，篩選得出最佳模型及預測結(jié)果。

3 仿真驗證和結(jié)果分析

3.1 研究區(qū)概況及數(shù)據(jù)來源

陜西省商洛市山陽縣的中村鎮(zhèn)，因其地處秦嶺山下，山脈溝壑眾多，屬于中、低山地形，山體土石量多達180多萬方，占地高達80%以上，位于地勢差異較大的峽谷地區(qū)，地形地質(zhì)復雜，山體石量多，更易引發(fā)災害。同時也屬于長江流域漢江水系，地區(qū)水源較多，河流較多、尤其在夏秋季降雨量也較多，年平均降水量達到671～865毫米，如果連續(xù)降雨量大容易導致土質(zhì)疏松[17]，從而增加了地質(zhì)災害的安全隱患點。

參照《T/CAGHP 006-2018泥石流災害防治工程勘察規(guī)范》、《滑坡崩塌泥石流災害調(diào)查規(guī)范(1：50 000)(DZ/T0261-2014)》[18]結(jié)合山陽縣實地監(jiān)測區(qū)域監(jiān)測泥石流活動數(shù)據(jù)，對泥石流發(fā)育機制及成災特征分析，本研究區(qū)域按照水源和物源成因劃分為崩塌型泥石流，其中固體物質(zhì)主要由滑坡崩塌等重力侵蝕提供[19-20]。去掉規(guī)范量級評分表中5分以下影響因子，最后結(jié)合監(jiān)測區(qū)實際泥石流數(shù)據(jù)得出11個影響因子，分別為溝岸山坡坡度、降雨量(24 h、1 h、10 min最大降雨量)、土壤含水率、溝床平均坡度、巖性影響、流域相對高差、河溝堵塞程度、河溝縱坡、產(chǎn)沙區(qū)溝槽橫斷面、松散物平均厚度、流域面積、泥砂沿程補給長度比、孔隙水壓力、沿溝松散物量、區(qū)域構(gòu)造影響、流域植被覆蓋率。因各降雨量24 h、1 h、10 min最大降雨量對泥石流的發(fā)生有極大的影響，所以選取暴雨強度R作為泥石流災害模型的影響因子。暴雨強度R計算如式(27)，各參數(shù)選取如表1所示。

表1 泥石流降雨量因子參數(shù)

R=K(H24/H24(D)+H1/H1(FD)+H1/6/H1/6(D))

(27)

3.2 數(shù)據(jù)預處理

監(jiān)測數(shù)據(jù)由于環(huán)境的影響會出現(xiàn)一些如缺失、離群或維度不統(tǒng)一的數(shù)據(jù)，這些數(shù)據(jù)對于模型的建立有極大的消極影響，會產(chǎn)生跳躍，且無法與其他數(shù)據(jù)統(tǒng)一，因此需要對監(jiān)測數(shù)據(jù)進行預處理。

1)異常值處理：監(jiān)測數(shù)據(jù)中存在一部分偏離傳感器本身范圍的值或偏離觀測值較大的值，不處理會影響數(shù)據(jù)預測準確性，距離達到5倍或者與均值的距離≥3倍標準差的數(shù)據(jù)稱之為離群點。

2)缺失值的處理：監(jiān)測數(shù)據(jù)通過泥石流災害區(qū)域多傳感器實時傳輸，傳輸過程中經(jīng)常會出現(xiàn)遺漏或者個別離群點情況，會出現(xiàn)失真損失有效信息，導致屬性值缺失不準確。按照屬性因素方法進行統(tǒng)計得出缺失率，本文劃分兩種類別數(shù)據(jù)的缺失值，如表2所示。

表2 數(shù)據(jù)缺失值

3)數(shù)據(jù)歸一化：監(jiān)測數(shù)據(jù)種類較多樣且數(shù)量較多，多傳感器數(shù)據(jù)量綱不同有較大的差異，使用原始數(shù)據(jù)直接建模對于預測的準確性有極大的影響，所以需要對數(shù)據(jù)進行歸一化處理，歸一化處理公式如(28)：

(28)

式(28)中，R為某因素歸一化處理后的數(shù)據(jù)，Rmin和Rmax表示某因素數(shù)據(jù)中的最小值及最大值。

3.3 數(shù)據(jù)影響因子篩選

由于樣本影響因子彼此之間存在相關性，為避免相關性對預測結(jié)果的影響，本文通過KPCA核主成分分析法選取成災因子，各主成分的特征值及貢獻率如圖4所示，實驗表明前6個主成分的累計的貢獻率已經(jīng)達到95.48%，覆蓋的信息超過了90%，覆蓋率達到要求，所以文中選取前6個影響因子作為泥石流災害模型訓練的輸入數(shù)據(jù)。并依據(jù)《T/CAGHP 006-2018泥石流災害防治工程勘察規(guī)范》及泥石流相關資料分析，得出影響因子與泥石流發(fā)生量化等級關系如表3所示。以陜西省山陽縣重點地災監(jiān)測區(qū)的歷史數(shù)據(jù)作為研究樣本，數(shù)據(jù)使用山陽縣2019年4月到2021年4月的10個監(jiān)測點的數(shù)據(jù)作為數(shù)據(jù)集，經(jīng)過數(shù)據(jù)預處理及成災因子選取后數(shù)據(jù)集總共篩選出1 300組數(shù)據(jù)，分別分為80%測試集和20%的兩個驗證集。

表3 泥石流影響因子量化表

圖4 影響因子特征值及主成分貢獻率

注:A.泥石流發(fā)生嚴重;B.泥石流中等發(fā)生;C.泥石流輕微發(fā)生;D.泥石流不發(fā)生。

3.4 超參數(shù)尋優(yōu)

LSSVM建模過程中調(diào)優(yōu)參數(shù)為正則化系數(shù)和核函數(shù)φ參數(shù)，文中選取SSA尋優(yōu)算法與遺傳算法 (GA，genetic algorithm)及網(wǎng)格搜索(GC，gridsearchCV)在相同1 040組訓練集對LSSVM模型的正則化系數(shù)γ和核函數(shù)φ參數(shù)進行尋優(yōu)。種群的規(guī)模設置為90，最大迭代次數(shù)設置為200，每個優(yōu)化算法分別進行60次獨立實驗，并分別畫出最優(yōu)適應度函數(shù)值與迭代次數(shù)曲線圖進行比對，結(jié)果如圖5所示，適應度函數(shù)值隨著迭代次數(shù)的增加而逐漸減小，最終搜尋到最優(yōu)參數(shù)后收斂。GC在第16次迭代大幅下降。跳出了局部最優(yōu)狀態(tài)，GA整個迭代過程收斂速度較慢，但也逐漸趨向最優(yōu)，SSA優(yōu)化效果最好，明顯引導種群向最優(yōu)位置處，說明使用SSA算法尋優(yōu)，對松鼠移動的位置不斷調(diào)整可以跳出局部最優(yōu)值，且收斂速度快，而且早熟現(xiàn)象明顯，能夠取得更小的適應度。最終選取正則化系數(shù)γ=0.274和核函數(shù)φ=7.642。

圖5 尋優(yōu)適應度曲線對比圖

3.5 仿真驗證及結(jié)果分析

為驗證模型的預測精度，引入模型評價指標AUC值，其為ROC曲線結(jié)合坐標軸圍成的面積值，范圍一般介于[0.5～1]之間，預測的真實性取決與AUC值接近1的程度，靠近1真實性高反之則反。MAE預測值真實誤差，RMSE預測值和真實值間偏離程度，MSE真實值與預測值差異，越接近零預測精度越高，如式(29)所示：

(29)

為驗證本研究優(yōu)化模型的準確性，將經(jīng)過數(shù)據(jù)預處理及降維后的訓練數(shù)據(jù)作為泥石流預測模型構(gòu)建的輸入數(shù)據(jù)，總共1 040組訓練集構(gòu)建泥石流預測模型，并通過10%驗證集1驗證各模型的準確性。實驗采用LSSVM作為泥石流災害預測模型，并用SSA算法超參數(shù)尋優(yōu)。使用同一個驗證集驗證未優(yōu)化的LSSVM模型及其他尋優(yōu)算法對LSSVM預測效果比對。利用預測結(jié)果計算模型的MAE、MSE和RMSE，值越接近于零精度越高，可以看出SSA-LSSVM的MAE、MSE和RMSE最小且接近于零，對比評估指標結(jié)果如圖6所示，傳統(tǒng)的LSSVM相對誤差較大，最大相對誤差達到1.72%，而SSA-LSSVM最大誤差達到0.19%，誤差是最低的，進一步說明了該模型預測的精度較高。

圖6 模型預測評估指標

為進一步驗證模型的穩(wěn)健性，選取10%驗證集2，將其打亂的130個監(jiān)測數(shù)據(jù)作為模型預測概率及預測等級誤差的評價，圖7為各模型尋優(yōu)LSSVM模型后的實際發(fā)生概率與預測發(fā)生概率比對圖。SSA尋優(yōu)后實際值與預測值基本吻合，擬合情況較好，極限的幾個數(shù)據(jù)27、38、89及111發(fā)生概率存在一些差異，但是其對應的風險預報等級與實際數(shù)據(jù)風險等級結(jié)果吻合，不影響預報的等級，多個算法模型在同一預測集上的預測等級結(jié)果如圖8所示，按照泥石流發(fā)生等級準確率降序排列：SSA達到100%，GA-LSSVM達到92.3%，GC-LSSVM達到90%。實驗說明引入SSA對LSSVM參數(shù)尋優(yōu)，泥石流發(fā)生的概率及等級預測準確率皆有明顯的提升。使用尋優(yōu)后的SSA-LSSVM模型對研究區(qū)域泥石流進行預測，從發(fā)生的概率及預測的等級兩方面都證明該模型具有較高的可行性。

圖7 模型預測預測比對圖

圖8 發(fā)生概率等級預測圖

此外通過AUC公式計算得出20%驗證集中各模型的ROC曲線如圖9所示，ROC曲線中橫坐標為假陽性率(FPR/1-Specificity特異度)，縱坐標為真陽性率(TPR/Sensitivity)，可以根據(jù)ROC曲線的面積下的AUC值看出各個預測模型對應評價指標的好壞，AUC值越高說明模型精度越佳，各模型AUC值均高于0.88，但是SSA-LSSVM模型指標更加，無論從測試時間、AUC值及ROC曲線均明顯優(yōu)于GA和GC尋優(yōu)模型，模型AUC均值為0.932，預測效果較其他模型理想。各模型的對比指標參數(shù)如表4所示，SSA-LSSVM模型相比其它模型平均測試事件最短且平均AUC值最高且接近于1。多組實驗數(shù)據(jù)結(jié)果證明SSA-LSSVM模型具有較好的預測效果，在泥石流災害預測中有較好的預測能力。

表4 模型參數(shù)及結(jié)果比對

圖9 模型ROC曲線

4 結(jié)束語

本文以普適的山陽縣中村鎮(zhèn)區(qū)域泥石流為例，結(jié)合泥石流全域的地形地貌成災機理，提出KPCA-SSA-LSSVM泥石流預測模型，在研究區(qū)實踐應用效果良好，說明模型具有一定的可行性和有效性。因此，可以得出以下結(jié)論：

1)參照《T/CAGHP 006-2018泥石流災害防治工程勘察規(guī)范》、《滑坡崩塌泥石流災害調(diào)查規(guī)范(1：50 000)(DZ/T0261-2014)》并結(jié)合山陽縣實地監(jiān)測區(qū)域，監(jiān)測泥石流活動數(shù)據(jù)，對泥石流發(fā)育機制及成災特征分析，選出11個成災因子，并使用KPCA主成分分析法依據(jù)因子的貢獻率篩選出6個重要的成災因子；

2)LSSVM建模過程中調(diào)優(yōu)參數(shù)為正則化系數(shù)和核函數(shù)參數(shù)，選取SSA尋優(yōu)算法與遺傳算法 (GA，genetic algorithm)及網(wǎng)格搜索(GC，GridSearchCV)在相同1 040組訓練集對LSSVM模型參數(shù)進行尋優(yōu)，解決參數(shù)隨機導致的精度不佳問題及陷入局部最優(yōu)問題。

3)將尋優(yōu)后的SSA-LSSVM預測結(jié)果與GA、GC參數(shù)尋優(yōu)模型預測結(jié)果比對，從AUC值、MAE、MSE、RMSE評價指標都驗證了SSA-LSSVM預測的精度。

4)使用尋優(yōu)后的SSA-LSSVM模型對研究區(qū)域泥石流進行預測，從發(fā)生的概率及預測的等級兩方面都證明該模型具有較高的可行性。