高澤民 丁明濤 楊國輝 黃 濤 張曉宇 周云濤 席傳杰
(①西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院, 成都 611756, 中國)
(②中鐵第一勘察設(shè)計院集團有限公司, 西安 710043, 中國)
受地震、降雨和人類活動等誘發(fā)因素的影響,泥石流災(zāi)害的暴發(fā)具有運動過程劇烈和破壞性強的特點,使其成為普遍而又危險的災(zāi)害類型之一。對于世界各地已經(jīng)發(fā)生的泥石流事件,在災(zāi)后常常造成嚴重的生命傷亡和經(jīng)濟損失(Hungr et al.,1987; Jakob et al.,2005; Nadim et al.,2006; 游勇等, 2011; Perov et al.,2017; Klime? et al.,2020; 劉傳正等, 2020)。特別是在地質(zhì)構(gòu)造、地形條件和氣候變化極其復(fù)雜的青藏高原及周緣地區(qū),以泥石流為代表的典型地質(zhì)災(zāi)害事件尤為突出,不僅發(fā)生的頻率高,而且破壞規(guī)模也較為顯著(陳寧生等, 2006; 崔鵬, 2009; 崔鵬等, 2010; 唐川等, 2011; Ni et al.,2012; Qiu et al.,2017; 魏昌利等, 2019)。而穿行于青藏高原之上的川藏鐵路正處于上述地質(zhì)環(huán)境中,工程建設(shè)在災(zāi)害脅迫下受到嚴峻的挑戰(zhàn)。因此,在空間尺度上合理量化和預(yù)測泥石流的潛在危險性和危險區(qū)域?qū)τ谖磥泶ú罔F路災(zāi)害防治、降低人類工程活動與自然環(huán)境競爭中的損失具有積極的作用和意義。
地質(zhì)災(zāi)害的危險性是指某種災(zāi)害在一定時期、一定的區(qū)域內(nèi)發(fā)生一定破壞規(guī)模事件的可能性,這種可能性可用空間上的概率來定量評估。更確切的說,危險性不僅要考慮災(zāi)害發(fā)生破壞的時間和頻率,也要考慮其預(yù)期的發(fā)生規(guī)模,是包含未來情況的一種評估(Brabb, 1984; Guzzetti, 2006)。為了平衡人-地競爭中存在的矛盾,最大限度地降低災(zāi)害損失,許多研究深入并有效地開展了防災(zāi)、減災(zāi)方面的工作。特別是在地質(zhì)災(zāi)害評估方面,建立了評估模型和指標(biāo)體系,提出了相關(guān)的評估技術(shù)手段,形成了比較廣泛和系統(tǒng)的災(zāi)害評估理論和方法。這些研究所涵蓋的區(qū)域覆蓋世界各地,用于評估的災(zāi)害類型豐富,取得了顯著的成果。主要的研究領(lǐng)域涉及全球或中小區(qū)域尺度的災(zāi)害風(fēng)險評價(劉艷輝等, 2019; Zou et al.,2019; Thiery et al.,2020)、易損性評價(Ding et al.,2016; 劉艷輝等, 2018; Nahayo et al., 2019; Gao et al.,2021)、易發(fā)性評價(楊光等, 2019; Orhan et al.,2020; Rodrigues et al.,2021)和危險性評價(Deyle et al.,1998; Guzzetti, 2006; 朱崇浩等, 2020)。所涉及的模型包括確定性、統(tǒng)計性和人工智能算法(Tan et al.,2015; Ding et al.,2016; Chen et al.,2017; Ahmed et al.,2020)。將這些方法互相結(jié)合并利用GIS技術(shù),即可初步形成地質(zhì)災(zāi)害評價技術(shù)體系中的相應(yīng)評估模塊。在上述模型中,基于貝葉斯優(yōu)化算法的隨機森林(The Bayesian optimization algorithm based random forest, TBOR)和梯度提升樹模型(The Bayesian optimization algorithm based gradient boosting tree, TBOG)分別為機器學(xué)習(xí)中兩種重要的集成算法,目前已經(jīng)被廣泛、成功地應(yīng)用并且顯示出了良好的災(zāi)害預(yù)測能力(Wang et al.,2015; Hong et al.,2016; Sahin, 2020; 劉艷輝等, 2021)。然而,對于泥石流溝分布較為廣泛的川藏鐵路孜熱—波密段,仍然未能將集成的機器學(xué)習(xí)算法應(yīng)用于該區(qū)域的災(zāi)害危險性評估,在區(qū)域災(zāi)害評估工作方法上存在一定的滯后和欠缺。
基于上述研究背景,為了完善和補充川藏鐵路孜熱—波密段地質(zhì)災(zāi)害危險性評估工作的內(nèi)容和方法,研究利用基于貝葉斯優(yōu)化算法的隨機森林和梯度提升樹兩種模型,選取區(qū)域典型的泥石流為災(zāi)害類型,結(jié)合11個地質(zhì)環(huán)境因子——高程、地表起伏度、坡度、坡向、平面曲率、剖面曲率、地表覆蓋類型、工程巖組、距河流距離、距斷層距離、距鐵路線路距離來預(yù)測泥石流危險性水平。該研究成果可以為川藏鐵路孜熱—波密線路防災(zāi)減災(zāi)工程建設(shè)提供參考和指導(dǎo)。
川藏鐵路孜熱—波密段地處青藏高原東緣,行政范圍隸屬西藏自治區(qū)林芝市波密縣境內(nèi)(圖1)。線路沿帕隆藏布江而下,東起林芝市波密縣,跨越通麥天險,西至林芝市巴宜區(qū)孜熱,全長約200km。區(qū)內(nèi)山地聚落人口分布分散,有國家干線拉林公路(國道318線)橫貫工作區(qū),交通設(shè)施相對落后。
圖1 川藏鐵路孜熱—波密段區(qū)域位置及泥石流災(zāi)害發(fā)育概況
川藏鐵路孜熱—波密段穿越念青唐古拉山東段、雅魯藏布江斷裂與三江斷裂應(yīng)力緩沖地帶,多期次構(gòu)造變形運動使巖體擠壓破碎產(chǎn)生大量節(jié)理和斷層等結(jié)構(gòu)面。其中, 對孜熱—波密段線路影響較大的為嘉黎斷裂和雅魯藏布江結(jié)合帶南、北界斷裂。區(qū)內(nèi)地形地貌復(fù)雜多樣,地勢起伏顯著,局部懸崖陡立,崎嶇難行。地勢總體上呈北高南低、東高西低態(tài)勢。主要地層為岡底斯-騰沖地層,出露基巖以念青唐古拉群片麻巖為主,局部有花崗巖、花崗閃長巖等侵入巖,分布在帕龍藏布江兩岸高山區(qū)域。該區(qū)地殼上升較為強烈,巖石受強風(fēng)化作用破碎明顯,河流侵蝕快速,夷平面破碎,多形成扇狀礫石堆積,易發(fā)生崩塌、滑坡、泥石流等地質(zhì)災(zāi)害(李佶航等, 2020)。
根據(jù)現(xiàn)場調(diào)查,研究區(qū)各泥石流溝分布于帕隆藏布流域河谷兩側(cè),溝道形態(tài)特征主要以中小流域、窄陡溝道型為主,溝口高程約為2600~2700m,主溝長度多介于4.88~7.17km,溝道兩側(cè)地形坡度多處于25°~55°,平均縱坡降多介于300‰~519‰之間,溝域面積范圍5.23~20.79km2,溝域河谷相對高差范圍約2500~2800m,具備泥石流暴發(fā)的地形條件。溝域內(nèi)海拔3800m以上區(qū)域分布有季節(jié)性冰雪區(qū),凍融侵蝕作用強烈,以下則以林地、灌木和草地為主。受印度洋西南季風(fēng)影響,孜熱—波密段所處的藏東南地區(qū)降水量較為充沛,具有春夏末降水量大、秋冬變小的特點。據(jù)波密站氣象資料記載,泥石流溝口位置海拔一帶的年平均氣溫約為8.8℃,年均降雨量約為898.6mm,高山區(qū)一帶平均氣溫1.0℃,年降雨量約為2400mm。氣候帶的垂直分帶差異使得溝內(nèi)冰雪積累、消融和地表徑流的匯集來源更加廣泛。結(jié)合現(xiàn)場勘察情況顯示,區(qū)內(nèi)泥石流溝的水源條件主要以冰雪融水、降雨和冰雪融水-降水融合型為主。此外,由于該區(qū)海洋性冰川活動強烈,流域內(nèi)溝道物源除滑坡、崩塌以及凍融風(fēng)化的溝道堆積物外,冰磧物、冰水堆積物也是重要的物源,堆積位置多處于海拔3700m以上的溝道中。冰磧物主要由礫石組成,大小混雜,呈次圓及次棱角狀,基質(zhì)為含礫亞黏土。溝道堆積物中的松散細粒物質(zhì)受到冰雪融水、降雨、洪水的裹挾和沖蝕作用,常使堆積塊石土具有一定的架空結(jié)構(gòu)特征,在不同的沖淤條件下影響著泥石流的堆積形態(tài),在出溝口形成結(jié)構(gòu)復(fù)雜的扇裝堆積特征。
隨機森林是機器學(xué)習(xí)中的一種集成算法之一,最早由Leo Breiman和Adele Cutler提出,用以解決數(shù)據(jù)的分類和回歸問題(Breiman, 2001)。它由多個決策樹構(gòu)建,每個決策樹通過不斷分裂將訓(xùn)練數(shù)據(jù)集分裂為兩個獨立的子數(shù)據(jù)集,最終形成樹形結(jié)構(gòu)。將完全樹中每一棵分類樹的預(yù)測結(jié)果結(jié)合就可得到隨機森林的預(yù)測結(jié)果。由于隨機森林具有高訓(xùn)練速度、通用性和極好的準確率等優(yōu)點,可以解決各類分類問題,因而在地質(zhì)災(zāi)害評價方面獲得較多成功的應(yīng)用(Zhang et al.,2017; 朱清華, 2019; 肖婷, 2020)。
梯度提升樹是機器學(xué)習(xí)中一種重要的集成算法。根據(jù)樣本的輸出特征值是否連續(xù),可將GBT模型分為分類和回歸兩類。在地質(zhì)災(zāi)害危險性評價中,由于模型的輸出特征值是離散的,因此對區(qū)域潛在的災(zāi)害發(fā)生概率預(yù)測和分級實質(zhì)上是分類問題。區(qū)別于回歸算法的是,分類算法無法直接從輸出類別中擬合類別輸出的誤差。為了解決這一問題,可以采用對數(shù)似然損失函數(shù)的方法。也就是說,可以用類別的概率預(yù)測值和真實概率值的差來擬合損失。對數(shù)似然損失函數(shù)又可分為二元分類和多元分類。在本文中,GBT分類算法僅僅涉及二元分類的對數(shù)似然損失函數(shù)。該模型可以表示為(Friedman,2002; Chen et al.,2016):
(1)
式中:M為迭代次數(shù),T(x,θm)為每次迭代產(chǎn)生的弱分類器,θm為損失函數(shù),可由式(2)計算:
(2)
式中:Fm-1(xi)為當(dāng)前迭代。
機器學(xué)習(xí)模型的參數(shù)優(yōu)化可以提升應(yīng)用模型的分類預(yù)測精度。貝葉斯優(yōu)化的目的在于求某個目標(biāo)函數(shù)的最小值,然后對目標(biāo)函數(shù)循環(huán)計算的評價結(jié)果建立另一個相應(yīng)的替代函數(shù)來尋求最小化目標(biāo)函數(shù)的值(Rong et al.,2020)。相比于傳統(tǒng)的網(wǎng)格搜索法需要手動調(diào)參、隨機網(wǎng)格搜索時間長的缺點,貝葉斯超參數(shù)優(yōu)化算法可以自動調(diào)整超參數(shù)搜索,并在每一次迭代過程中重新選擇參數(shù)時參考循環(huán)之前的評估結(jié)果,較大幅度縮短了尋優(yōu)時間,有效提升優(yōu)化效率。優(yōu)化后的TBOR和TBOG模型的超參數(shù)尋優(yōu)結(jié)果可見表1。
表1 基于貝葉斯優(yōu)化的TBOR和TBOG危險性評價模型超參數(shù)調(diào)整結(jié)果
泥石流災(zāi)害的發(fā)生與地球內(nèi)動力作用和外在環(huán)境條件密切相關(guān)。本研究共選擇了11個危險性評價因子作為模型輸入特征,所涵蓋的范圍涉及地質(zhì)、地理和環(huán)境條件等方面,能客觀地反映并提供可靠的危險性基礎(chǔ)信息,并通過對輸入特征參數(shù)閾值的進一步劃分,得到了泥石流災(zāi)害發(fā)育的分區(qū)統(tǒng)計特征(表2)。
表2 TBOR和TBOG危險性評價模型輸入特征分區(qū)數(shù)據(jù)統(tǒng)計
高程:從災(zāi)害數(shù)量來看,高程在1596~3125m范圍內(nèi)的泥石流事件最多,有80處,占災(zāi)害總量的46.512%。其次是區(qū)間3126~3743m,發(fā)育有35處泥石流災(zāi)害,占20.349%。881~7176m區(qū)間內(nèi)僅發(fā)育有7次事件,占泥石流災(zāi)害總量的4.07%(表2)。總體來看,隨著高程的升高,地質(zhì)災(zāi)害發(fā)生數(shù)量整體呈現(xiàn)下降的趨勢。
地表起伏:該因子表示研究區(qū)內(nèi)最大高程和最小高程的差值,反映地表起伏的陡緩。本研究利用的地表起伏值由空間分辨為120m×120m的柵格單元求高差值得到。然后按照自然斷點法將高差重分類為5級區(qū)間。在這些區(qū)間范圍中,災(zāi)害主要發(fā)生于0~81m分組內(nèi),共發(fā)生52次事件,發(fā)生頻率30.233%。217~303m和304~1910m區(qū)間內(nèi)災(zāi)害發(fā)生頻率較低,頻率比分別僅為13.372%和2.907%。區(qū)間82~303m內(nèi)發(fā)生泥石流事件115次,發(fā)生頻率居中,占比為66.9%(表2)。
坡度:研究區(qū)地形坡度主要分布于0°~81°范圍內(nèi)。其中: 0°~11°區(qū)間內(nèi)災(zāi)害數(shù)量最多,為47處, 40°~81°范圍內(nèi)發(fā)育最少的11處災(zāi)害。22°~30°和32°~39°兩個區(qū)間內(nèi)發(fā)育的泥石流數(shù)量接近,分別為39次和40次(表2)。
坡向:研究區(qū)災(zāi)害點坡向主要分布于216°~287°和144°~215°范圍之內(nèi),分別發(fā)生50次和40次災(zāi)害事件。其余區(qū)間內(nèi)災(zāi)害點坡向分布范圍介于18°~34°之間(表2)。
平面曲率:在該輸入特征的5級區(qū)間中, 216°~287°和144°~215°范圍內(nèi)的災(zāi)害點數(shù)量最多,共有88個,占總災(zāi)害點數(shù)量的51.16%。其余3個區(qū)間災(zāi)害發(fā)生的頻率介于13.953%~20.930%(表2)。
剖面曲率:將坡度值再進行一次坡度計算可得到剖面曲率值,其值域為0~15.1。其中:災(zāi)害數(shù)發(fā)育最多的區(qū)間為1.3~2.4,而5.8~15.1則最少,僅發(fā)生有1次。0~1.2、2.5~3.9和4~5.7區(qū)間災(zāi)害發(fā)生頻率居中(表2)。
土地利用類型:如表2所示,研究區(qū)地質(zhì)災(zāi)害主要集中于高山林地和草地,分別有74處和46處,占災(zāi)害總數(shù)的43.023%和26.744%。其次為冰川和永久積雪,發(fā)生比例為17.442%。而其余地表覆蓋類型的地質(zhì)災(zāi)害發(fā)生數(shù)量僅有22處,占總災(zāi)害數(shù)量的12.791%。
工程巖組:不用時期不同類型的巖石抗侵蝕能力有較大差異,對災(zāi)害發(fā)育影響顯著。區(qū)內(nèi)地質(zhì)災(zāi)害主要分布于變質(zhì)巖組、碎屑巖組和變質(zhì)巖組內(nèi),災(zāi)害數(shù)量分別為74處、44處和32處,三者所占的比例分別達43.023%、36.661%和12.209%(表2)。而其余各個地質(zhì)時期內(nèi)災(zāi)害發(fā)生的數(shù)量相對較少。
距線路距離:該特征參數(shù)用于衡量災(zāi)害點距鐵路線路的距離的影響。從研究建立的5級緩沖區(qū)統(tǒng)計數(shù)據(jù)可以看出,區(qū)間0~200m內(nèi)的災(zāi)害發(fā)育數(shù)量和頻率分別為158個和91.860%(表2)。其他4個區(qū)間內(nèi)地質(zhì)災(zāi)害數(shù)量分布較少。
距斷層距離:利用緩沖分析功能可得到地質(zhì)災(zāi)害分布與距斷層距離間的空間關(guān)系(表2)。統(tǒng)計結(jié)果顯示,地質(zhì)災(zāi)害發(fā)生總量和頻率在0~200m的區(qū)間內(nèi)最高,占整個研究區(qū)的92.442%,其他區(qū)間內(nèi)地質(zhì)災(zāi)害數(shù)量分布不明顯。
距河流距離:該特征參數(shù)的5級緩沖區(qū)由河流分布矢量數(shù)據(jù)獲得。由表2可知,地質(zhì)災(zāi)害主要分布在距河流0~200m的范圍內(nèi),共123處,占地質(zhì)災(zāi)害總量的71.512%。200~400m地質(zhì)災(zāi)害數(shù)量為38處,其他3個區(qū)間內(nèi)的地質(zhì)災(zāi)害總數(shù)量為11處??傮w而言,泥石流災(zāi)害數(shù)量和災(zāi)害密度均具有隨著距河流距離的增加而減小的趨勢。
研究基于Python計算平臺實現(xiàn)基于貝葉斯優(yōu)化算法的隨機森林和梯度提升樹模型的構(gòu)建。首先對11個指標(biāo)因子的柵格及矢量值進行空間數(shù)據(jù)的提取。由于研究區(qū)范圍較大,并且受到計算性能的限制,我們將原始空間分辨率為30m×30m的柵格影像重采樣為120m×120m,得到744476個網(wǎng)格單元。這樣不僅使預(yù)測范圍單元化,而且提升運算效率,并且地形精度仍然相適應(yīng),可以較好地保證預(yù)測準確性。然后,利用GIS平臺將這些空間數(shù)據(jù)分布值提取到2970個網(wǎng)格中,得到每個評價單元所對應(yīng)的11個特征值和1個標(biāo)簽值。通過在Python中讀取上述采集到的樣本數(shù)據(jù)庫,構(gòu)建了維數(shù)為2970×11的矩陣。緊接著,按照7︰3的比例將這些樣本進行隨機分割,獲得矩陣維數(shù)為2079×11的訓(xùn)練集和891×11測試集數(shù)據(jù)。最后,分別構(gòu)建基于貝葉斯優(yōu)化算法的隨機森林和梯度提升樹分類器,完成對744476個網(wǎng)格單元泥石流危險性概率的計算。為了更直觀地表現(xiàn)出泥石流危險性的空間分布特征,研究將每個單元的概率預(yù)測結(jié)果賦值到每一個劃分的評價單元中,并采用自然斷點劃分方法,將危險性值分為高危險、較高危險、中等危險、較低危險和低危險5級區(qū)間,最終得到了2個評價模型的危險性區(qū)劃圖(圖2)。
圖2 由TBOR(a)和TBOG(b)模型計算得到的泥石流危險性分布圖
表3 TBOR和TBOG模型中不同危險性分區(qū)數(shù)據(jù)與地質(zhì)災(zāi)害分布特征統(tǒng)計
研究結(jié)果顯示(表3),在TBOR分類模型中,高危險區(qū)所占的面積比例最大,為38.862%,共記錄有104次泥石流事件,占區(qū)內(nèi)災(zāi)害總量的60.465%,地質(zhì)災(zāi)害密度僅為12.577處/(102km2)。較高-低危險性區(qū)間的面積比相對較低,分別為7.713%、21.655%、14.193%和17.577%,分別占0.581%、16.860%、9.302%和13.372%的災(zāi)害點密度,且不同分區(qū)內(nèi)的危險性空間分布特征多呈條帶狀。TBOG模型中各危險性分區(qū)所占的比例和TBOR模型的預(yù)測結(jié)果大致接近,面積比例最大的仍然為高危險性區(qū)間,占比為56.806%,災(zāi)害點密度12.940處/(102km2), 而中等危險區(qū)間比例最小,為7.828%。較高、較低和低危險性區(qū)間的占比則分別為9.774%、9.487%和16.105%,對應(yīng)16.279%、1.163%和12.209%的密度比例。綜合以上結(jié)果可以看出,地質(zhì)災(zāi)害高-較高危險區(qū)在深切河谷內(nèi)分布最為廣泛,沿水系及道路線狀地物展布明顯。其次,受復(fù)雜的構(gòu)造地貌格局和氣候條件等孕災(zāi)環(huán)境的影響和控制,加之工程活動對地質(zhì)條件的擾動,泥石流物源、溝道特征等要素可能在部分線路段存在相關(guān)的改造。因此,該線路段安全建設(shè)仍然處于一定的災(zāi)害脅迫下,應(yīng)特別重視災(zāi)害防治和應(yīng)急保障工作。
預(yù)測模型綜合分類的效果可由受試者工作特征曲線(receiver operating characteristic curve, ROC)來檢驗。為了對分類效果給一個確定的評價值,ROC曲線下方與坐標(biāo)軸圍成的面積(area under curve, AUC)值被廣泛采用(Beguería, 2006; Vakhshoori et al.,2018)。在基于貝葉斯優(yōu)化的隨機森林和梯度提升樹測試集中,將模型各自的預(yù)測標(biāo)簽與實際標(biāo)簽值進行對比,就可以得到真陽率(True positive rate, TPR)和假陽率(False positive rate, TPR)這兩種分類預(yù)測比率。并分別以TPR為縱軸,F(xiàn)PR為橫軸繪制每一個樣本點的預(yù)測比率,即可得出ROC曲線。同時,ROC曲線上(0, 0)和(1, 1)兩點間的連線代表分類閾值(Chance)。通常情況下,如果ROC曲線在該閾值線的上方,并且曲線越靠近(0, 1)這個點,那么曲線下方包含的面積值(AUC)越大,表征模型分類效果越好。當(dāng)AUC值大于0.9時,則模型具有極高的預(yù)測精度(Fawcett, 2006)。根據(jù)本研究計算結(jié)果(圖3),TBOR和TBOG模型分類器的AUC值分別為0.89和0.83,即GPM模型的預(yù)測精度分別為89%和83%,顯示出兩種模型在用于泥石流危險性預(yù)測時具有較高的可靠性。同時,TBOR模型的ROC曲線位置明顯更高,表明了該算法具有更高的分類精度。
圖3 TBOR和TBOG模型的ROC-AUC對比結(jié)果
研究在闡述地質(zhì)災(zāi)害危險性評價研究現(xiàn)狀的基礎(chǔ)上,通過對實際勘察資料數(shù)據(jù)的整理分析,運用基于貝葉斯優(yōu)化算法的隨機森林和梯度提升樹模型,建立了川藏鐵路孜熱—波密段泥石流災(zāi)害評價模型及指標(biāo)體系,完成了該區(qū)的泥石流災(zāi)害危險性區(qū)劃判定。主要得到以下結(jié)論:
(1)依據(jù)TBOR和TBOG模型,結(jié)合研究區(qū)11個地質(zhì)環(huán)境特征指標(biāo),得到了研究區(qū)泥石流災(zāi)害危險性分區(qū)圖。其中, 高-較高危險區(qū)的比例分別占研究區(qū)總面積的56.439%和66.580%,并呈現(xiàn)出帶狀分布特征。
(2)將泥石流災(zāi)害點與危險性區(qū)劃圖疊加后的結(jié)果顯示,TBOR和TBOG兩種模型在高-較高危險區(qū)內(nèi)泥石流災(zāi)害分布更為集中,各自累計共發(fā)生有127和128處災(zāi)害事件,災(zāi)害密度分別為12.577處/(102km2)和12.940處/(102km2),災(zāi)害密度比例分別達73.837%和74.419%。
(3)受限于收集的數(shù)據(jù)類型和精度差異以及評價模型的選擇,本論文得出的評價結(jié)果在線路部分區(qū)域仍然存在一定的不確定性。因此,該研究還需要結(jié)合野外實地調(diào)查情況,才能為線路防災(zāi)減災(zāi)防護工程建設(shè)提供更加科學(xué)的參考和指導(dǎo)。