川藏鐵路孜熱—波密段泥石流災害危險性評價*

高澤民 丁明濤 楊國輝 黃 濤 張曉宇 周云濤 席傳杰

(①西南交通大學地球科學與環境工程學院， 成都 611756， 中國)

(②中鐵第一勘察設計院集團有限公司， 西安 710043， 中國)

0 引 言

受地震、降雨和人類活動等誘發因素的影響，泥石流災害的暴發具有運動過程劇烈和破壞性強的特點，使其成為普遍而又危險的災害類型之一。對于世界各地已經發生的泥石流事件，在災后常常造成嚴重的生命傷亡和經濟損失(Hungr et al.，1987; Jakob et al.，2005; Nadim et al.，2006; 游勇等， 2011; Perov et al.，2017; Klime? et al.，2020; 劉傳正等， 2020)。特別是在地質構造、地形條件和氣候變化極其復雜的青藏高原及周緣地區，以泥石流為代表的典型地質災害事件尤為突出，不僅發生的頻率高，而且破壞規模也較為顯著(陳寧生等， 2006; 崔鵬， 2009; 崔鵬等， 2010; 唐川等， 2011; Ni et al.，2012; Qiu et al.，2017; 魏昌利等， 2019)。而穿行于青藏高原之上的川藏鐵路正處于上述地質環境中，工程建設在災害脅迫下受到嚴峻的挑戰。因此，在空間尺度上合理量化和預測泥石流的潛在危險性和危險區域對于未來川藏鐵路災害防治、降低人類工程活動與自然環境競爭中的損失具有積極的作用和意義。

地質災害的危險性是指某種災害在一定時期、一定的區域內發生一定破壞規模事件的可能性，這種可能性可用空間上的概率來定量評估。更確切的說，危險性不僅要考慮災害發生破壞的時間和頻率，也要考慮其預期的發生規模，是包含未來情況的一種評估(Brabb， 1984; Guzzetti， 2006)。為了平衡人-地競爭中存在的矛盾，最大限度地降低災害損失，許多研究深入并有效地開展了防災、減災方面的工作。特別是在地質災害評估方面，建立了評估模型和指標體系，提出了相關的評估技術手段，形成了比較廣泛和系統的災害評估理論和方法。這些研究所涵蓋的區域覆蓋世界各地，用于評估的災害類型豐富，取得了顯著的成果。主要的研究領域涉及全球或中小區域尺度的災害風險評價(劉艷輝等， 2019; Zou et al.，2019; Thiery et al.，2020)、易損性評價(Ding et al.，2016; 劉艷輝等， 2018; Nahayo et al.， 2019; Gao et al.，2021)、易發性評價(楊光等， 2019; Orhan et al.，2020; Rodrigues et al.，2021)和危險性評價(Deyle et al.，1998; Guzzetti， 2006; 朱崇浩等， 2020)。所涉及的模型包括確定性、統計性和人工智能算法(Tan et al.，2015; Ding et al.，2016; Chen et al.，2017; Ahmed et al.，2020)。將這些方法互相結合并利用GIS技術，即可初步形成地質災害評價技術體系中的相應評估模塊。在上述模型中，基于貝葉斯優化算法的隨機森林(The Bayesian optimization algorithm based random forest， TBOR)和梯度提升樹模型(The Bayesian optimization algorithm based gradient boosting tree， TBOG)分別為機器學習中兩種重要的集成算法，目前已經被廣泛、成功地應用并且顯示出了良好的災害預測能力(Wang et al.，2015; Hong et al.，2016; Sahin， 2020; 劉艷輝等， 2021)。然而，對于泥石流溝分布較為廣泛的川藏鐵路孜熱—波密段，仍然未能將集成的機器學習算法應用于該區域的災害危險性評估，在區域災害評估工作方法上存在一定的滯后和欠缺。

基于上述研究背景，為了完善和補充川藏鐵路孜熱—波密段地質災害危險性評估工作的內容和方法，研究利用基于貝葉斯優化算法的隨機森林和梯度提升樹兩種模型，選取區域典型的泥石流為災害類型，結合11個地質環境因子——高程、地表起伏度、坡度、坡向、平面曲率、剖面曲率、地表覆蓋類型、工程巖組、距河流距離、距斷層距離、距鐵路線路距離來預測泥石流危險性水平。該研究成果可以為川藏鐵路孜熱—波密線路防災減災工程建設提供參考和指導。

1 區域概況

1.1 地理位置

川藏鐵路孜熱—波密段地處青藏高原東緣，行政范圍隸屬西藏自治區林芝市波密縣境內(圖1)。線路沿帕隆藏布江而下，東起林芝市波密縣，跨越通麥天險，西至林芝市巴宜區孜熱，全長約200km。區內山地聚落人口分布分散，有國家干線拉林公路(國道318線)橫貫工作區，交通設施相對落后。

圖1 川藏鐵路孜熱—波密段區域位置及泥石流災害發育概況

1.2 地質地形地貌

川藏鐵路孜熱—波密段穿越念青唐古拉山東段、雅魯藏布江斷裂與三江斷裂應力緩沖地帶，多期次構造變形運動使巖體擠壓破碎產生大量節理和斷層等結構面。其中， 對孜熱—波密段線路影響較大的為嘉黎斷裂和雅魯藏布江結合帶南、北界斷裂。區內地形地貌復雜多樣，地勢起伏顯著，局部懸崖陡立，崎嶇難行。地勢總體上呈北高南低、東高西低態勢。主要地層為岡底斯-騰沖地層，出露基巖以念青唐古拉群片麻巖為主，局部有花崗巖、花崗閃長巖等侵入巖，分布在帕龍藏布江兩岸高山區域。該區地殼上升較為強烈，巖石受強風化作用破碎明顯，河流侵蝕快速，夷平面破碎，多形成扇狀礫石堆積，易發生崩塌、滑坡、泥石流等地質災害(李佶航等， 2020)。

1.3 泥石流形成條件及基本特征分析

根據現場調查，研究區各泥石流溝分布于帕隆藏布流域河谷兩側，溝道形態特征主要以中小流域、窄陡溝道型為主，溝口高程約為2600～2700m，主溝長度多介于4.88～7.17km，溝道兩側地形坡度多處于25°～55°，平均縱坡降多介于300‰～519‰之間，溝域面積范圍5.23～20.79km2，溝域河谷相對高差范圍約2500～2800m，具備泥石流暴發的地形條件。溝域內海拔3800m以上區域分布有季節性冰雪區，凍融侵蝕作用強烈，以下則以林地、灌木和草地為主。受印度洋西南季風影響，孜熱—波密段所處的藏東南地區降水量較為充沛，具有春夏末降水量大、秋冬變小的特點。據波密站氣象資料記載，泥石流溝口位置海拔一帶的年平均氣溫約為8.8℃，年均降雨量約為898.6mm，高山區一帶平均氣溫1.0℃，年降雨量約為2400mm。氣候帶的垂直分帶差異使得溝內冰雪積累、消融和地表徑流的匯集來源更加廣泛。結合現場勘察情況顯示，區內泥石流溝的水源條件主要以冰雪融水、降雨和冰雪融水-降水融合型為主。此外，由于該區海洋性冰川活動強烈，流域內溝道物源除滑坡、崩塌以及凍融風化的溝道堆積物外，冰磧物、冰水堆積物也是重要的物源，堆積位置多處于海拔3700m以上的溝道中。冰磧物主要由礫石組成，大小混雜，呈次圓及次棱角狀，基質為含礫亞黏土。溝道堆積物中的松散細粒物質受到冰雪融水、降雨、洪水的裹挾和沖蝕作用，常使堆積塊石土具有一定的架空結構特征，在不同的沖淤條件下影響著泥石流的堆積形態，在出溝口形成結構復雜的扇裝堆積特征。

2 研究方法

2.1 隨機森林模型

隨機森林是機器學習中的一種集成算法之一，最早由Leo Breiman和Adele Cutler提出，用以解決數據的分類和回歸問題(Breiman， 2001)。它由多個決策樹構建，每個決策樹通過不斷分裂將訓練數據集分裂為兩個獨立的子數據集，最終形成樹形結構。將完全樹中每一棵分類樹的預測結果結合就可得到隨機森林的預測結果。由于隨機森林具有高訓練速度、通用性和極好的準確率等優點，可以解決各類分類問題，因而在地質災害評價方面獲得較多成功的應用(Zhang et al.，2017; 朱清華， 2019; 肖婷， 2020)。

2.2 梯度提升樹模型

梯度提升樹是機器學習中一種重要的集成算法。根據樣本的輸出特征值是否連續，可將GBT模型分為分類和回歸兩類。在地質災害危險性評價中，由于模型的輸出特征值是離散的，因此對區域潛在的災害發生概率預測和分級實質上是分類問題。區別于回歸算法的是，分類算法無法直接從輸出類別中擬合類別輸出的誤差。為了解決這一問題，可以采用對數似然損失函數的方法。也就是說，可以用類別的概率預測值和真實概率值的差來擬合損失。對數似然損失函數又可分為二元分類和多元分類。在本文中，GBT分類算法僅僅涉及二元分類的對數似然損失函數。該模型可以表示為(Friedman，2002; Chen et al.，2016)：

(1)

式中：M為迭代次數，T(x，θm)為每次迭代產生的弱分類器，θm為損失函數，可由式(2)計算：

(2)

式中：Fm-1(xi)為當前迭代。

2.3 貝葉斯優化

機器學習模型的參數優化可以提升應用模型的分類預測精度。貝葉斯優化的目的在于求某個目標函數的最小值，然后對目標函數循環計算的評價結果建立另一個相應的替代函數來尋求最小化目標函數的值(Rong et al.，2020)。相比于傳統的網格搜索法需要手動調參、隨機網格搜索時間長的缺點，貝葉斯超參數優化算法可以自動調整超參數搜索，并在每一次迭代過程中重新選擇參數時參考循環之前的評估結果，較大幅度縮短了尋優時間，有效提升優化效率。優化后的TBOR和TBOG模型的超參數尋優結果可見表1。

表1 基于貝葉斯優化的TBOR和TBOG危險性評價模型超參數調整結果

3 模型輸入特征空間發育規律

泥石流災害的發生與地球內動力作用和外在環境條件密切相關。本研究共選擇了11個危險性評價因子作為模型輸入特征，所涵蓋的范圍涉及地質、地理和環境條件等方面，能客觀地反映并提供可靠的危險性基礎信息，并通過對輸入特征參數閾值的進一步劃分，得到了泥石流災害發育的分區統計特征(表2)。

表2 TBOR和TBOG危險性評價模型輸入特征分區數據統計

高程：從災害數量來看，高程在1596～3125m范圍內的泥石流事件最多，有80處，占災害總量的46.512%。其次是區間3126～3743m，發育有35處泥石流災害，占20.349%。881～7176m區間內僅發育有7次事件，占泥石流災害總量的4.07%(表2)?？傮w來看，隨著高程的升高，地質災害發生數量整體呈現下降的趨勢。

地表起伏：該因子表示研究區內最大高程和最小高程的差值，反映地表起伏的陡緩。本研究利用的地表起伏值由空間分辨為120m×120m的柵格單元求高差值得到。然后按照自然斷點法將高差重分類為5級區間。在這些區間范圍中，災害主要發生于0～81m分組內，共發生52次事件，發生頻率30.233%。217～303m和304～1910m區間內災害發生頻率較低，頻率比分別僅為13.372%和2.907%。區間82～303m內發生泥石流事件115次，發生頻率居中，占比為66.9%(表2)。

坡度：研究區地形坡度主要分布于0°～81°范圍內。其中： 0°～11°區間內災害數量最多，為47處， 40°～81°范圍內發育最少的11處災害。22°～30°和32°～39°兩個區間內發育的泥石流數量接近，分別為39次和40次(表2)。

坡向：研究區災害點坡向主要分布于216°～287°和144°～215°范圍之內，分別發生50次和40次災害事件。其余區間內災害點坡向分布范圍介于18°～34°之間(表2)。

平面曲率：在該輸入特征的5級區間中， 216°～287°和144°～215°范圍內的災害點數量最多，共有88個，占總災害點數量的51.16%。其余3個區間災害發生的頻率介于13.953%～20.930%(表2)。

剖面曲率：將坡度值再進行一次坡度計算可得到剖面曲率值，其值域為0～15.1。其中：災害數發育最多的區間為1.3～2.4，而5.8～15.1則最少，僅發生有1次。0～1.2、2.5～3.9和4～5.7區間災害發生頻率居中(表2)。

土地利用類型：如表2所示，研究區地質災害主要集中于高山林地和草地，分別有74處和46處，占災害總數的43.023%和26.744%。其次為冰川和永久積雪，發生比例為17.442%。而其余地表覆蓋類型的地質災害發生數量僅有22處，占總災害數量的12.791%。

工程巖組：不用時期不同類型的巖石抗侵蝕能力有較大差異，對災害發育影響顯著。區內地質災害主要分布于變質巖組、碎屑巖組和變質巖組內，災害數量分別為74處、44處和32處，三者所占的比例分別達43.023%、36.661%和12.209%(表2)。而其余各個地質時期內災害發生的數量相對較少。

距線路距離：該特征參數用于衡量災害點距鐵路線路的距離的影響。從研究建立的5級緩沖區統計數據可以看出，區間0～200m內的災害發育數量和頻率分別為158個和91.860%(表2)。其他4個區間內地質災害數量分布較少。

距斷層距離：利用緩沖分析功能可得到地質災害分布與距斷層距離間的空間關系(表2)。統計結果顯示，地質災害發生總量和頻率在0～200m的區間內最高，占整個研究區的92.442%，其他區間內地質災害數量分布不明顯。

距河流距離：該特征參數的5級緩沖區由河流分布矢量數據獲得。由表2可知，地質災害主要分布在距河流0～200m的范圍內，共123處，占地質災害總量的71.512%。200～400m地質災害數量為38處，其他3個區間內的地質災害總數量為11處?？傮w而言，泥石流災害數量和災害密度均具有隨著距河流距離的增加而減小的趨勢。

4 結果分析與驗證

研究基于Python計算平臺實現基于貝葉斯優化算法的隨機森林和梯度提升樹模型的構建。首先對11個指標因子的柵格及矢量值進行空間數據的提取。由于研究區范圍較大，并且受到計算性能的限制，我們將原始空間分辨率為30m×30m的柵格影像重采樣為120m×120m，得到744476個網格單元。這樣不僅使預測范圍單元化，而且提升運算效率，并且地形精度仍然相適應，可以較好地保證預測準確性。然后，利用GIS平臺將這些空間數據分布值提取到2970個網格中，得到每個評價單元所對應的11個特征值和1個標簽值。通過在Python中讀取上述采集到的樣本數據庫，構建了維數為2970×11的矩陣。緊接著，按照7︰3的比例將這些樣本進行隨機分割，獲得矩陣維數為2079×11的訓練集和891×11測試集數據。最后，分別構建基于貝葉斯優化算法的隨機森林和梯度提升樹分類器，完成對744476個網格單元泥石流危險性概率的計算。為了更直觀地表現出泥石流危險性的空間分布特征，研究將每個單元的概率預測結果賦值到每一個劃分的評價單元中，并采用自然斷點劃分方法，將危險性值分為高危險、較高危險、中等危險、較低危險和低危險5級區間，最終得到了2個評價模型的危險性區劃圖(圖2)。

圖2 由TBOR(a)和TBOG(b)模型計算得到的泥石流危險性分布圖

表3 TBOR和TBOG模型中不同危險性分區數據與地質災害分布特征統計

研究結果顯示(表3)，在TBOR分類模型中，高危險區所占的面積比例最大，為38.862%，共記錄有104次泥石流事件，占區內災害總量的60.465%，地質災害密度僅為12.577處/(102km2)。較高-低危險性區間的面積比相對較低，分別為7.713%、21.655%、14.193%和17.577%，分別占0.581%、16.860%、9.302%和13.372%的災害點密度，且不同分區內的危險性空間分布特征多呈條帶狀。TBOG模型中各危險性分區所占的比例和TBOR模型的預測結果大致接近，面積比例最大的仍然為高危險性區間，占比為56.806%，災害點密度12.940處/(102km2)， 而中等危險區間比例最小，為7.828%。較高、較低和低危險性區間的占比則分別為9.774%、9.487%和16.105%，對應16.279%、1.163%和12.209%的密度比例。綜合以上結果可以看出，地質災害高-較高危險區在深切河谷內分布最為廣泛，沿水系及道路線狀地物展布明顯。其次，受復雜的構造地貌格局和氣候條件等孕災環境的影響和控制，加之工程活動對地質條件的擾動，泥石流物源、溝道特征等要素可能在部分線路段存在相關的改造。因此，該線路段安全建設仍然處于一定的災害脅迫下，應特別重視災害防治和應急保障工作。

預測模型綜合分類的效果可由受試者工作特征曲線(receiver operating characteristic curve， ROC)來檢驗。為了對分類效果給一個確定的評價值，ROC曲線下方與坐標軸圍成的面積(area under curve， AUC)值被廣泛采用(Beguería， 2006; Vakhshoori et al.，2018)。在基于貝葉斯優化的隨機森林和梯度提升樹測試集中，將模型各自的預測標簽與實際標簽值進行對比，就可以得到真陽率(True positive rate， TPR)和假陽率(False positive rate， TPR)這兩種分類預測比率。并分別以TPR為縱軸，FPR為橫軸繪制每一個樣本點的預測比率，即可得出ROC曲線。同時，ROC曲線上(0， 0)和(1， 1)兩點間的連線代表分類閾值(Chance)。通常情況下，如果ROC曲線在該閾值線的上方，并且曲線越靠近(0， 1)這個點，那么曲線下方包含的面積值(AUC)越大，表征模型分類效果越好。當AUC值大于0.9時，則模型具有極高的預測精度(Fawcett， 2006)。根據本研究計算結果(圖3)，TBOR和TBOG模型分類器的AUC值分別為0.89和0.83，即GPM模型的預測精度分別為89%和83%，顯示出兩種模型在用于泥石流危險性預測時具有較高的可靠性。同時，TBOR模型的ROC曲線位置明顯更高，表明了該算法具有更高的分類精度。

圖3 TBOR和TBOG模型的ROC-AUC對比結果

5 結 論

研究在闡述地質災害危險性評價研究現狀的基礎上，通過對實際勘察資料數據的整理分析，運用基于貝葉斯優化算法的隨機森林和梯度提升樹模型，建立了川藏鐵路孜熱—波密段泥石流災害評價模型及指標體系，完成了該區的泥石流災害危險性區劃判定。主要得到以下結論：

(1)依據TBOR和TBOG模型，結合研究區11個地質環境特征指標，得到了研究區泥石流災害危險性分區圖。其中， 高-較高危險區的比例分別占研究區總面積的56.439%和66.580%，并呈現出帶狀分布特征。

(2)將泥石流災害點與危險性區劃圖疊加后的結果顯示，TBOR和TBOG兩種模型在高-較高危險區內泥石流災害分布更為集中，各自累計共發生有127和128處災害事件，災害密度分別為12.577處/(102km2)和12.940處/(102km2)，災害密度比例分別達73.837%和74.419%。

(3)受限于收集的數據類型和精度差異以及評價模型的選擇，本論文得出的評價結果在線路部分區域仍然存在一定的不確定性。因此，該研究還需要結合野外實地調查情況，才能為線路防災減災防護工程建設提供更加科學的參考和指導。