尚艷芳， 李麗敏， 溫宗周， 王朝陽， 夏夢凡

(西安工程大學(xué)電子信息學(xué)院， 西安 710600)

中國地形地貌較為復(fù)雜，其中山區(qū)地形所占比例占70%，且處于地震多發(fā)帶，再加之不合理的工程活動，使中國山區(qū)居民飽受泥石流災(zāi)害的困擾，泥石流災(zāi)害破壞力巨大，容易導(dǎo)致大量人員傷亡和經(jīng)濟(jì)損失。因此，采用科學(xué)、先進(jìn)的泥石流危險(xiǎn)性評價方法，建立有效的預(yù)測預(yù)報(bào)系統(tǒng)是非常必要的。

針對泥石流危險(xiǎn)性評價的研究開始較早，已取得較為豐富的成果，并且，隨著泥石流的研究不斷深入，各種其他學(xué)科理論和方法不斷引入泥石流研究領(lǐng)域，中外許多學(xué)者對泥石流危險(xiǎn)性評價也有了更透徹認(rèn)識。李麗敏等[1]針對泥石流多因素影響，提出最小二乘支持向量機(jī)的泥石流災(zāi)害預(yù)報(bào)模型，但模型輸入多維數(shù)據(jù)，計(jì)算機(jī)占用較多時間資源和內(nèi)存資源，造成網(wǎng)絡(luò)迭代速度較慢。根據(jù)泥石流形成條件不同，鄧恩松等[2]針對兩大不同類型的(降雨型和冰川型)泥石流進(jìn)行研究，以中巴公路奧布段為研究區(qū)，結(jié)合了主成分分析法和灰色系統(tǒng)理論建立評估模型，但此種研究方法只適用于降雨型和冰川型泥石流災(zāi)害預(yù)測，具有一定的局限性。王俊豪等[3]使用層次分析法來確定泥石流危險(xiǎn)性評價中各影響因素的權(quán)重，通過引用多位專家的評判標(biāo)準(zhǔn)減小主觀誤差對結(jié)果的影響，但沒有進(jìn)行進(jìn)一步的驗(yàn)證性分析，可能導(dǎo)致多變量之間相互疊加，影響預(yù)報(bào)準(zhǔn)確性。王艷錦等[4]使用特征金字塔網(wǎng)絡(luò)對成康鐵路沿線泥石流進(jìn)行危險(xiǎn)性評估，研究結(jié)果表明了機(jī)器學(xué)習(xí)理論在突發(fā)性地質(zhì)災(zāi)害領(lǐng)域的有效性和可行性。張曉東[5]以寧夏鹽池縣作為地質(zhì)災(zāi)害研究區(qū)，使用遙感和多源信息融合技術(shù)進(jìn)行泥石流危險(xiǎn)性評價，這一方法雖然效果較好，但由于需要借助遙感技術(shù)，因此不具有普適性。王晨輝等[6]采用支持向量機(jī)模型預(yù)測泥石流的危險(xiǎn)程度，有大量數(shù)據(jù)可以供模型訓(xùn)練學(xué)習(xí)，這種方法雖然精度得到了保證但模型建立過程復(fù)雜，計(jì)算量較大。

針對上述預(yù)測方法中存在的問題，現(xiàn)采用主成分分析(principal component analysis，PCA)線性降維算法，對數(shù)據(jù)集做特征壓縮；之后，使用核極限學(xué)習(xí)機(jī)(kernel based extreme learning machine，KELM) 預(yù)測泥石流危險(xiǎn)性，并利用螢火蟲算法(firefly algorithm，F(xiàn)A)對其進(jìn)行參數(shù)尋優(yōu)，從而得到泥石流危險(xiǎn)性的分類結(jié)果；最后，將本文方法與核極限學(xué)習(xí)機(jī)模型、多分類支持向量機(jī)模型和BP(error back propagation，BP)神經(jīng)網(wǎng)絡(luò)模型進(jìn)行比較，驗(yàn)證FA-KELM模型的優(yōu)越性。

1 研究方法介紹

1.1 主成分分析算法(PCA)

主成分分析(PCA)采用幾個互相正交的變量信息代替多源的信息量，目的是以最少的新變量揭示原可觀測指標(biāo)，在數(shù)據(jù)特征壓縮的同時，最大程度避免特征信息間重疊，盡量減少原指標(biāo)包含信息的損失[7-8]。包含如下步驟。

(1)設(shè)i組j維變量構(gòu)成的數(shù)據(jù)矩陣為

(1)

式(1)中：xij為j維變量的第i個樣本，將xij標(biāo)準(zhǔn)化處理，得

(2)

(2)計(jì)算相關(guān)系數(shù)矩陣R。

R=[rij]m×n

(3)

(3)計(jì)算累計(jì)貢獻(xiàn)率βi。

(4)

式(4)中：λj為特征矩陣R對應(yīng)的特征值。

(4)選定主成分個數(shù)。一般情況下，當(dāng)累計(jì)貢獻(xiàn)率不少于85%時，則選定有m維主成分[9]，與之對應(yīng)的有m個特征值Λm=diag[λ1，λ2，…，λm]和m個特征向量Wm=[w1，w2，…，wm]。

1.2 核極限學(xué)習(xí)機(jī)(KELM)

核極限學(xué)習(xí)機(jī)(KELM)和其他傳統(tǒng)前饋神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)類似，但在訓(xùn)練學(xué)習(xí)時，所需迭代參數(shù)少，學(xué)習(xí)速度快，其輸出樣本函數(shù)F(x)可用矩陣[10]表示為

F(x)=h(x)β=Hβ=L

(5)

β=H*L

(6)

式中:x為給定的模型輸入因子；h(x)、H為隱含層神經(jīng)元的輸出；β為隱含層和輸出層之間的連接權(quán)重；L為網(wǎng)絡(luò)學(xué)習(xí)后的期望輸出；H*為H的廣義逆矩陣。通過求解最小二乘解[式(7)]可以得到輸出權(quán)重矩陣的結(jié)果。

(7)

式(7)中：I為單位矩陣；C為正則化系數(shù)；本文研究中的核函數(shù)選用高斯核函數(shù)，即

(8)

式(8)中：g為核參數(shù)；xi、xj為模型輸入因子，引入核函數(shù)后，核矩陣為

QELM=h(xi)·h(xj)=K(xi，xj)

(9)

則表達(dá)式(5)網(wǎng)絡(luò)輸出可描述為

(10)

1.3 基于FA算法訓(xùn)練KELM模型

雖然KELM算法克服了傳統(tǒng)ELM算法的缺點(diǎn)，但基于核的ELM則存在模型參數(shù)選擇問題，本文研究中通過螢火蟲算法[11-12](firefly algorithm，F(xiàn)A)調(diào)整KELM模型中的超參數(shù)，解決因人工設(shè)置參數(shù)而導(dǎo)致預(yù)測結(jié)果的隨機(jī)性和不確定性，提高建模精度和預(yù)測性能。該模型的訓(xùn)練過程如表1所示。

(11)

(12)

rij=‖xi-xj‖

(13)

式中:I0為最大亮度；m0為最大吸引度；γ為光強(qiáng)吸收系數(shù)；xi和xj分別為螢火蟲個體i和j的所在位置。個體i向另一個個體j空間轉(zhuǎn)移的依據(jù)公式為

表1 FA算法訓(xùn)練KELM模型過程Table 1 FA algorithm training KELM model process

(14)

2 研究區(qū)概況

北川縣坐落于四川省盆地西北部向川西高原的過渡帶上，峰巒起伏，山地自然景觀占全縣一半以上的面積，溝壑縱橫，季風(fēng)氣候明顯，陰雨天氣較多，年均降雨量1 399.1 mm，空間分布不均，具有東南向西北變小的規(guī)律[13]，同時，區(qū)內(nèi)最高海拔4 769 m，最低點(diǎn)540 m，不穩(wěn)定溝床比有明顯規(guī)律，易于泥石流輸移，流域切割密度較大，使得固體物源更易于參與泥石流活動。根據(jù)唐川等[14]、汪月鵑[15]調(diào)查，該區(qū)域內(nèi)相對高差超過1 000 m，溝谷斜坡坡度多為15°～40°，地形險(xiǎn)峻，地層巖性多樣，復(fù)雜的地質(zhì)構(gòu)造為流域的支溝發(fā)育和泥石流的發(fā)生提供有利條件(圖1)。

3 建模過程和結(jié)果分析

3.1 數(shù)據(jù)預(yù)處理

本文中數(shù)據(jù)摘自四川省北川縣國家重點(diǎn)地質(zhì)災(zāi)害監(jiān)測項(xiàng)目的歷史數(shù)據(jù)，選取其中72條泥石流溝基礎(chǔ)數(shù)據(jù)作為研究樣本。綜合《泥石流危險(xiǎn)性評價》[16]等，最終選擇8種初始影響因子，分別為流域切割密度X1、流域面積X2、流域相對高差X3、不穩(wěn)定溝床比X4、松散固體物質(zhì)儲量X5、年均降雨量X6、地震烈度X7和主溝長度X8。

為了克服輸入影響因子丟失或者異常值對地質(zhì)災(zāi)害預(yù)報(bào)準(zhǔn)確性的影響，分別針對噪聲大、存在野值、數(shù)據(jù)種類多的問題，使用以下3種方法進(jìn)行解決。

(1)適當(dāng)濾波法。本文使用移動平均濾波器對原始數(shù)據(jù)進(jìn)行平滑處理，使用smooth函數(shù)，設(shè)置移動平均濾波器的默認(rèn)窗寬為5。

圖1 北川縣泥石流分布圖Fig.1 Debris flow distribution map in Beichuan County

(15)

若Fji>F(n，a)，則判定該數(shù)據(jù)屬于異常值，應(yīng)該丟棄對應(yīng)的數(shù)據(jù)xji。

(3)歸一化。在本數(shù)據(jù)集中有8種數(shù)據(jù)，為避免對后續(xù)變量的運(yùn)算和分析造成很大的干擾，加快模型的求解速度，得到合理的結(jié)果，將所有綱量映射至[0，1]區(qū)間。轉(zhuǎn)換函數(shù)公式為

(16)

式(16)中：x′為經(jīng)過歸一化后的數(shù)據(jù)；x為樣本原始數(shù)據(jù)；xmax、xmin為樣本數(shù)據(jù)中的最大值、最小值。

為了驗(yàn)證預(yù)處理方法在本模型中的有效性，將預(yù)處理前后的數(shù)據(jù)分別作為FA-KELM預(yù)測模型的輸入，對比了數(shù)據(jù)預(yù)處理前后的分類結(jié)果，如圖2 所示。

從圖2對比結(jié)果可知，通過必要的預(yù)處理工作，剔除了異常值，數(shù)據(jù)的取值范圍均在客觀合理范圍內(nèi)，模型預(yù)測準(zhǔn)確率達(dá)到90.47%；然而，把不規(guī)范的數(shù)據(jù)作為輸入時，其預(yù)測準(zhǔn)確率僅有71.42%。仿真結(jié)果表明，在數(shù)據(jù)輸入模型前，確保數(shù)據(jù)格式一致化以及數(shù)據(jù)的完整性是十分有必要的。

圖2 數(shù)據(jù)預(yù)處理前后FA-KELM模型預(yù)測結(jié)果對比Fig.2 Comparison of FA-KELM model prediction results before and after data pretreatment

3.2 PCA降維

在數(shù)據(jù)集中隨機(jī)選擇51條泥石流溝數(shù)據(jù)作為網(wǎng)絡(luò)模型訓(xùn)練，剩余數(shù)據(jù)用于驗(yàn)證模型。選擇累計(jì)貢獻(xiàn)率不少于85%的成分作為主成分，利用PCA算法選擇出4種綜合指標(biāo)，涵蓋了8個原始指標(biāo)體系的絕大部分信息，達(dá)到了數(shù)據(jù)信息精煉化，提取后的變量如圖3所示，將線性組合后的新變量輸入到FA-KELM模型中進(jìn)行訓(xùn)練學(xué)習(xí)并驗(yàn)證，降維前后的測試時間分別為100.69、82.47 ms。有效提高KELM模型的學(xué)習(xí)能力。

由圖3中PCA主成分提取變量結(jié)果可得，前4個成分的累計(jì)貢獻(xiàn)率為45.428%、63.297%、76.132%、86.259%，且到第4個成分時的累計(jì)貢獻(xiàn)率達(dá)到了86.259%，大于定義主成分臨界值85%。因此利用PCA主成分分析法能夠?qū)⒊跏?維變量數(shù)據(jù)降至4維，模型數(shù)據(jù)結(jié)構(gòu)復(fù)雜度控制到更低且消除了各因子之間的相關(guān)冗余度。所得到的主成分系數(shù)矩陣如表2所示。

根據(jù)表2中主成分系數(shù)矩陣可得，原始變量和主成分之間的線性關(guān)系表達(dá)式為

(17)

表2 主成分系數(shù)矩陣Table 2 Principal component coefficient matrix

圖3 PCA主成分提取變量Fig.3 PCA principal component extraction variables

式(17)中：F1～F4為降維后的4維主成分變量；X1～X8為原始8維影響因子數(shù)據(jù)。

3.3 仿真驗(yàn)證及結(jié)果分析

3.3.1 參數(shù)選擇

KELM算法性能優(yōu)劣取決于懲罰系數(shù)C和核函數(shù)參數(shù)g的選擇，利用FA優(yōu)化算法對以上參數(shù)進(jìn)行優(yōu)化選擇，分類器的懲罰系數(shù)和核參數(shù)優(yōu)化范圍設(shè)置為(C，g)∈[2-5，210]，為避免實(shí)驗(yàn)的偶然性且計(jì)算量適中，種群規(guī)模n設(shè)置為35，最大迭代次數(shù)為50，得到的最優(yōu)參數(shù)C和g分別為142.413和1.802，優(yōu)化參數(shù)花費(fèi)時間為21.62 s，圖4為尋優(yōu)算法適應(yīng)度值曲線圖，螢火蟲算法以極快的收斂速度取得更高的適應(yīng)度函數(shù)值，尋優(yōu)效率高，整體收斂能力強(qiáng)?；赑ython3.7以及MATLAB2017b進(jìn)行編程實(shí)現(xiàn)網(wǎng)絡(luò)模型和優(yōu)化算法，CPU為Intel(R) Core(TM) i7-10700F CPU@2.90 GHz，RAM為16.0 GB。

圖4 尋優(yōu)算法適應(yīng)度曲線Fig.4 A chart of the adaptability curve of the optimization algorithm

3.3.2 評估模型訓(xùn)練

首先通過PCA對泥石流數(shù)據(jù)做特征提取，將訓(xùn)練集和測試集按3∶1劃分，并且將危險(xiǎn)性劃分為4個等級，分別為: 低度危險(xiǎn)、中度危險(xiǎn)、高度危險(xiǎn)、極高危險(xiǎn)，其級別分別對應(yīng)1級、2級、3級、4級?；贔A-KELM的泥石流危險(xiǎn)性預(yù)測模型的流程如圖5所示。初始時，使用隨機(jī)產(chǎn)生的參數(shù)作為FA的初始種群，將訓(xùn)練集輸入FA-KELM網(wǎng)絡(luò)進(jìn)行訓(xùn)練。預(yù)測準(zhǔn)確率作為個體適應(yīng)度函數(shù)，為去除較差子群，更新螢火蟲位置，以最大適應(yīng)度函數(shù)為標(biāo)準(zhǔn)選擇，使得模型獲得最佳調(diào)優(yōu)的參數(shù)，最后利用測試集測試該模型的準(zhǔn)確性。

3.3.3 性能評估

選取在泥石流預(yù)測方面應(yīng)用較為廣泛的核極限學(xué)習(xí)機(jī)[17](kernel based extreme learning machine，KELM)、多分類支持向量機(jī)[18-19](multi-class support vector machine， MSVM)以及誤差逆?zhèn)鞑ニ惴╗20-21](error back propagation，BP)作為對比模型，并采用混淆矩陣、曲線下面積(area under curve， AUC)、ROC曲線(receiver operating characteristic curve， ROC)評價分類器性能，在相同條件下FA-KELM模型與預(yù)測模型進(jìn)行效果對比。圖6和圖7給出了各對比模型在相同數(shù)據(jù)集下泥石流危險(xiǎn)性分類的混淆矩陣和ROC曲線。

圖5 基于FA-KELM的泥石流危險(xiǎn)性評價流程圖Fig.5 Flow chart of debris flow hazard assessment based on FA-KELM

由圖6對比可知，本文提出的模型在對角線上的值均高于其他模型在對角線上的值，得到最佳混淆矩陣，KELM預(yù)測性能稍好，MSVM算法的預(yù)測標(biāo)簽與真實(shí)標(biāo)簽之間存在較明顯的出入，BP算法預(yù)測結(jié)果偏離實(shí)際泥石流危險(xiǎn)性等級，并且預(yù)測模型的優(yōu)劣與ROC曲線下方面積的大小有關(guān)，面積越大，說明模型表現(xiàn)出色，根據(jù)圖7的ROC曲線結(jié)果可知，F(xiàn)A-KELM曲線的面積為最大，對角線代表分類能力為0的一條線，ROC曲線離對角線越遠(yuǎn)，AUC越大，表明模型的分類能力越強(qiáng)。各模型具體AUC如表3所示。

根據(jù)表3對比可知，本文提出的FA-KELM模型在ROC曲線下面積AUC均值為90.6%，因此，與其他網(wǎng)絡(luò)相比，其擁有最大AUC 均值，說明該算法更優(yōu)，在泥石流危險(xiǎn)性評估中具有較強(qiáng)競爭優(yōu)勢，有效提高分類效果。

圖6 不同模型結(jié)果的混淆矩陣對比圖Fig.6 Comparison of confusion matrix of the prediction results by different methods

圖7 不同模型ROC曲線對比Fig.7 Comparison of ROC curves of different models

表3 4種模型AUC均值對比Table 3 Comparison of AUC mean for 4 models

4 結(jié)論

以泥石流危害程度為研究對象，建立了基于PCA算法和經(jīng)過FA算法參數(shù)優(yōu)化的KELM地質(zhì)災(zāi)害評估模型。利用四川省北川縣監(jiān)測點(diǎn)數(shù)據(jù)進(jìn)行實(shí)驗(yàn)驗(yàn)證，同時，將本文方法與未進(jìn)行優(yōu)化改進(jìn)的KELM模型、多分類支持向量機(jī)模型、BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型的輸出結(jié)果進(jìn)行比較，通過對比驗(yàn)證，表明本文選用的FA-KELM模型在泥石流危險(xiǎn)性評估方面效果更優(yōu)。具體結(jié)果如下。

(1) 泥石流災(zāi)害的成因是復(fù)雜的，致災(zāi)因子較多，因此利用PCA進(jìn)行數(shù)據(jù)變換和降維處理，經(jīng)過特征選擇的數(shù)據(jù)更適應(yīng)模型的訓(xùn)練，去除了因子之間的冗余信息，降低了網(wǎng)絡(luò)運(yùn)算量，消耗更少運(yùn)算時間資源，提高模型準(zhǔn)確性。

(2) 由于數(shù)據(jù)是在復(fù)雜嚴(yán)峻的外部環(huán)境下采集到的，不可避免地存在不符合要求、重復(fù)、錯誤數(shù)據(jù)等現(xiàn)象，對輸出結(jié)果會產(chǎn)生影響，因此對數(shù)據(jù)集進(jìn)行了預(yù)處理，將預(yù)處理前后的數(shù)據(jù)分別作為模型的輸入，得到不同的結(jié)果，驗(yàn)證了預(yù)處理工作是模型訓(xùn)練前的重要環(huán)節(jié)，經(jīng)過數(shù)據(jù)處理后，提高了數(shù)據(jù)質(zhì)量，模型的準(zhǔn)確率得到了提高。

(3)將FA優(yōu)化的KELM方法應(yīng)用于泥石流危險(xiǎn)性評估中。選用FA優(yōu)化算法對KELM中的參數(shù)尋優(yōu)，并通過ROC曲線和對應(yīng)的AUC均值對模型進(jìn)行評估效果的評價，相較于其他3種傳統(tǒng)算法，F(xiàn)A-KELM算法展現(xiàn)出更好的效果，評估結(jié)果與實(shí)際情況基本一致，在實(shí)際工程中是可行有效的。

(4)地質(zhì)災(zāi)害的成功預(yù)報(bào)是極其復(fù)雜的系統(tǒng)工程，本文研究主要針對泥石流的危險(xiǎn)性進(jìn)行評估，而對于災(zāi)害發(fā)生強(qiáng)度、范圍和時間的預(yù)測也是工程中的難題。同時泥石流災(zāi)害從孕育到致災(zāi)過程受多種因素導(dǎo)致，且受地形地貌、氣象水文影響巨大，本文提出的模型提供一定的參考價值，但還需進(jìn)一步提高與完善，以便更好地運(yùn)用于其他流域地區(qū)的實(shí)際工程中。