林芝市公路地質災害易發性區劃及氣象預報模型

陶星宇, 柳錦寶, 黃志剛, 巴 桑, 卓 瑪, 陳 軍, 肖 斌, 高瑜蓮, 費曉燕

(1.成都信息工程大學 資源環境學院, 四川 成都 610225; 2.福建省氣象服務中心,福建 福州 350001; 3.西藏高原大氣環境科學研究所, 西藏 拉薩 850000; 4.西藏自治區氣候中心,西藏 拉薩 850000; 5.北京師范大學 地理學部, 北京 100875; 6.四川公眾項目咨詢管理有限公司, 四川 成都 610051)

林芝市境內主要公路大多分布于河流兩側的階地上，使得因強降水引起的地質災害的威脅以沖毀公路為主，其中，降水誘發的滑坡、崩塌、泥石流是林芝市地質災害中威脅、損失最大的災種[1]。截止到2019年，通過市縣地質災害調查與區劃查明的地質災害點約152處，隱患點約1 213處，而由降水誘發的地質災害占全市各類地質災害總數的80%，幾乎每年都會發生幾十次不同程度、不同規模的地質災害，造成重大損失。近年來，地質災害的發生愈加頻繁，造成的經濟損失越來越重，2016年10月11日，林芝市易貢茶場白龍溝發生大型泥石流地質災害，造成直接經濟損失2 154.59萬元；2018年10月17日，雅魯藏布江林芝市米林縣加拉村段發生滑坡地質災害，造成直接經濟損失118 158.94萬元。在汛期，林芝市大部分地質條件比較復雜的路段往往會發生滑坡、崩塌、泥石流等嚴重的地質災害，對沿線交通設施、群眾的生命安全構成嚴重的威脅。因此，加強林芝市公路沿線地質災害防治及氣象預報工作具有重要意義。

預報模型是成功開展地質災害氣象預報的基礎，近年來，地質災害氣象預報模型的研究得到了國內外學者及政府的高度重視[2-7]，劉傳正等[8]針對三峽庫區建立了降水判據模型，提出了地質災害預報與防治措施，給地質災害預報研究提供了新方向；Keiko[9]基于多元線性回歸法對誘發地質災害的降水量因子進行了分析[9]，這種只考慮降水量而未考慮到地質災害背景條件的閾值模型，結果差異很大，實際應用受到限制[10-14]。因此，對地質背景與降水量相結合的預報模型開始了研究，在多因子模型的基礎上，浙江、貴州、江蘇、廣東、四川、廣西、甘肅等省開始構建了地質災害氣象預報系統[15-21]；費曉燕、陳宮燕等[22-23]結合降水資料和地質環境特征，基于邏輯回歸方法分別對四川省和林芝市進行地質災害氣象預報研究。而現在國內大多研究的是大范圍地質災害氣象預報，小范圍預報和服務對象等方面考慮較少[24-25]，對于構建小范圍公路地質災害氣象預報模型方面的研究更是少之又少。目前，陳洪凱等[26]選取了暴雨、溫度、植被覆蓋率、災害體、年均降水量、巖性條件、路基位置、地震烈度、地貌類型、地質構造等10個因子，針對川藏公路提出了地質災害危險性區劃模型；孟兆興等[27]基于信息量法、邏輯回歸模型及空間邏輯回歸模型對川藏公路進行了易發性評價。目前，對于公路沿線地質災害的研究主要集中在易發性評價方面，未結合降水量因子進行分析，無法進行地質災害預報服務。因此，本研究則在地質災害易發性評價的基礎上，進一步進行地質災害氣象預報模型研究。本研究通過建立林芝市公路沿線15 km緩沖區，對林芝市公路沿線地質災害易發性區劃和降水誘發地質災害關系進行研究，構建林芝市公路沿線地質災害氣象預報模型，為實施林芝市境內公路地質災害氣象預報提供了依據。

1 研究區概況

本研究以林芝市公路沿線15 km緩沖區為研究區，林芝市位于西藏自治區的東南部喜馬拉雅山脈中段，南靠岡底斯山脈，北傍念青唐古拉山，地處雅魯藏布江中下游，屬高山峽谷區，山高谷深，土質基本為沙壤結構，森林覆蓋率為46%，是中國第三大林區。研究區內主要為山地地貌與河谷、山間盆地侵蝕堆積地貌區，地形完整性差，地形起伏大，平均海拔3 000 m，最低點海拔155 m，是陸地垂直地貌落差最大的地區。同時，該區域是強烈的構造隆起區與不同的地質構造區的交匯地帶[28]，該研究區屬于藏東南地區強烈的擠壓帶，包括喜馬拉雅板片、雅魯藏布江縫合帶、岡底斯—念青唐古拉板塊、班公錯—怒江縫合帶、羌塘—三江復合板塊、金沙江縫合帶[29-30]。研究區內均發育有各個時代的地層，巖石類型分布復雜，主要有泥盆系、三疊系的板巖、片巖、片麻巖、大理巖以及侏羅紀的泥巖等[31]。受印度洋暖流與北方寒流的影響，形成了熱帶、亞熱帶、溫帶和寒帶并存的復雜氣候帶，年降水量達到650 mm左右，降水多且集中。對研究區內2010—2019年日降水量數據進行年均月降水量統計分析，由于2015年氣象站點數據較少，不具備典型性，則將其剔除，剔除后發現5—9月(主汛期)的年均月降水量較大，其中，7月最大，年均月降水量達到了4 055 mm。因此，在地形地貌、地質構造、地層巖性和降水等因素的共同作用下使得林芝市成為西藏自治區地質災害最嚴重的地區，其中，2010—2019年研究區災害點及隱患點分布如圖1所示。

圖1 林芝市公路沿線災害點及隱患點分布

2 研究數據及方法

2.1 數據來源

本研究中使用的數據主要包括Landsat 8影像、DEM數據、土地利用、土壤、降水量、河流、道路、斷層、地質災害點及隱患點等數據。其中，Landsat 8影像源于美國地質勘探局官網USGS(https:∥earthexplorer.usgs.gov/)；DEM數據源于NASA官網(https:∥www.nasa.gov/)SRTM1 30 m數字高程數據；土地利用數據源于NASA官網2018年MOD12Q1產品數據；土壤類型數據源于《1∶100萬中華人民共和國土壤圖》；河流、道路數據源于Open Street Map官網(https://www.openstreetmap.org/)；斷層數據源于全國地質信息資料網(http://www.drc.cgs.gov.cn)全國1∶250萬地質數據；降水量、地質災害點及隱患點數據由林芝市氣象局及自然資源局提供，時間段為2010—2019年；坡度、溝谷密度數據是基于DEM數據在ArcGIS軟件中提取得到；裸巖率數據是基于Landsat8影像提取得到，時間為2019年11月29日。

2.2 研究方法

(1) 邏輯回歸方法。邏輯回歸方法在解決地質災害問題中，具有計算方法簡單、能自動篩選影響因子等優點。可將對地質災害發生產生影響的區劃因子作為自變量，將地質災害是否發生作為因變量(1代表發生，0代表未發生)[32]。考慮到各因子為非連續變量，本研究使用二元邏輯回歸方法進行分析[33]，計算公式為：

(1)

式中：P為地質災害發生的概率值； (1-P)為未發生地質災害的概率值；X1，X2，…，Xn為區劃因子；B1，B2，…，Bn為各區劃因子對應的邏輯回歸系數。

(2) 單因子信息量方法。由于各區劃因子有不同標準的數據類型和數據量級，不能直接對其進行分析。所以，在使用邏輯回歸方法之前需要將各區劃因子進行異質數據同化。信息量法在本質上是一種數據驅動的定量統計方法，2001年信息量法被用于地質災害易發性分級中。信息量法的主要思路為：根據歷史災害點，對各區劃因子進行分級，將各區劃因子的真實值轉化為反映該研究區區劃因子等級的信息量值，再根據信息量值來評判各區劃因子與地質災害的相關性。本研究即采用信息量法用于邏輯回歸模型的異質數據同化。

對于單一評價因素的信息量i，計算公式為：

(2)

式中：S為研究區內評價單元的總面積；N為評價單元內發生地質災害的總數量；Si為評價單元內特定類別的面積；Ni為分布在評價單元內特定類別的地質災害的數量。

3 地質災害易發性區劃研究

3.1 地質災害易發性區劃因子選取

通過前期對林芝市公路沿線15 km緩沖區內的地形地貌、地質環境和生態環境的研究，發現其與地質災害的分布緊密相關，同時，在中外地質災害易發性區劃模型的研究中，多數使用了坡度、斷層、降水量等因子作為區劃因子。因此，本研究最終選定8個區劃因子作為地質災害易發性區劃模型的評價指標，包括斷層密度、土地利用類型、土壤類型、坡度、溝谷密度、年均降水量、裸巖率和距河流的距離。再根據文獻調研、實際經驗和地質災害點及隱患點分布規律對8個區劃因子進行重分類(圖2)。

3.1.1 地形地貌指標

(1) 坡度。坡度與氣象水文條件、地表狀況等因素緊密相關，是評價地質災害的重要因素之一。一般情況下，坡度較大的地方會為地質災害的發生提供大量物質來源。但坡度增大到一定限度后，發生地質災害的概率將不會增大。本研究將坡度劃分為0°～10°，10°～20°，20°～30°，30°～40°，>40°共5級。

(2) 溝谷密度。溝谷密度是氣候、地形地貌、地層巖性、植被等因素綜合影響的反映，溝谷密度與地面破碎程度成正比，地面越破碎，地表物質的穩定性越低，地質災害的發生概率越高。本研究基于自然斷點法將溝谷密度分為0～0.12，0.12～0.22，0.22～0.32，>0.32 km/km2共4級。

(3) 裸巖率。裸巖率是劃分石漠化等級的重要參數，有些位居山地陡坡、土層較薄地帶極易發展為裸露或裸巖地，從而引發地質災害。本研究依據相等間隔法將裸巖率分為<0，0～0.2，0.2～0.4，0.4～0.6，0.6～0.8，0.8～1。

圖2 林芝市公路沿線各區劃因子分布

3.1.2 地質環境指標

(1) 斷層密度。斷層是影響地質災害發生的重要因素之一，總體上距離斷層越遠，地質災害點的分布密度越小，地質災害的強發育區是在斷層兩側的10 km范圍內[34]。一般采用距斷層的距離或斷層密度來判斷斷層對地質災害的影響程度。本研究選用斷層密度作為評價指標，基于自然斷點法將其分為0～0.015，0.015～0.05，0.05～0.085，0.085～0.12，>0.12 km/km2共5級。

(2) 距河流的距離。河流作為水文地質因素，其對地質巖層的影響十分顯著，由此引發的地質災害較為常見。水流反復流經會引起巖層塌陷，從而增加地質災害的發生概率。本研究根據研究區實際情況，基于ArcGIS多環緩沖區工具將距河流的距離分為0～1 000，1 000～2 000，2 000～3 000，3 000～4 000，>4 000 m共5級。

3.1.3 生態環境指標

(1) 土地利用類型。土地利用類型可體現地表植被覆蓋的狀況和人類工程活動的強度。其中，植被覆蓋狀況與地質災害發生密切相關，同時，人類工程活動(道路、房屋)也加劇了地質災害的發生。本研究基于MOD12Q1土地利用產品數據的類型說明表進行分類。

(2) 土壤類型。不同土壤的特性會對地質災害造成不同程度的影響，如土壤保水性能、滲水速率和含沙量等都會影響地質災害的發生。本研究基于1∶100萬中國土壤數據庫中的土類進行分類。

(3) 年均降水量。降水條件對地質災害的發生非常重要，降水下滲、洪水掏蝕河岸、水位變化等均會誘發地質災害，根據前期對林芝市公路沿線降水特征分析，絕大多數地質災害是由強降水誘發的。本研究基于自然斷點法將年均降水量分為462～716，716～834，834～966，966～1 163，>1 163 mm共5級。

為了避免各區劃因子數據類型和數據量級的不同而對結果產生影響，本研究基于單因子信息量方法對重分類后的8個指標進行異質數據同化，計算結果見表1。因子信息量值越大，表明地質災害易發性概率越大。

表1 單因子信息量計算結果

3.2 地質災害易發性區劃模型建立

本研究基于邏輯回歸方法建立地質災害易發性區劃模型，以1 213個有效隱患點作為建模數據，152個有效災害點作為模型驗證數據(圖1)。同時，為了使得二元邏輯回歸分析更具有科學性，本研究基于張錫濤、常鳴等人[35-36]的研究，在以隱患點為中心的3 km緩沖區外隨機選取對應數量的無災害發生樣本點。使用ArcGIS軟件提取隱患點樣本和無災害發生樣本的區劃因子信息量值。將地質災害是否發生作為因變量，斷層密度、土地利用、土壤類型、坡度、溝谷密度、年降水量、裸巖率和距河流的距離信息量值作為自變量，在SPSS軟件中構建二元邏輯回歸模型。為了減少無災害發生樣本數據對模型產生的影響，多次進行無災害發生樣本的選取并循環建模，最優結果的Nagelkerke方差可達到0.769，Cox & SnellR方差可達到0.577。表2為區劃因子邏輯回歸結果，其中，B為系數，S.E為標準差，Wald統計量越大表示該因子的重要性越強，從表2可以看出，各區劃因子重要性依次為：河流、坡度、土壤類型、斷層密度、年降水量、土地利用、溝谷密度、裸巖率，8個區劃因子的顯著性均小于0.05，顯著性較高，具有統計學意義。

表2 研究區區劃因子邏輯回歸方程

將表2中各因子系數代入公式中計算地質災害易發性概率值，計算公式為：

(3)

式中：H為地質災害易發性概率值；X1—X8分別為斷層密度、土壤類型、坡度、溝谷密度、年降水量、土地利用、裸巖率、距河流的距離信息量值。

將8個區劃因子的信息量值和系數代入公式(3)計算出地質災害易發性概率值，用ArcGIS軟件進行加權疊加分析得出林芝市公路地質災害易發性概率值，然后進行易發性分區，形成分區圖(圖3)。以0.1為間隔，將林芝市公路地質災害易發性概率分為10級。從圖3可知，林芝市公路地質災害易發性概率超過0.8的主體區域主要分布在朗縣中部、林芝縣東部、墨脫縣境內扎墨公路、察隅縣東部的察佐公路及察丙路附近；工布江縣、波密縣大部分區域易發性概率小于0.4，不易發生災害。總體趨勢是離道路越近，易發性概率越大。

圖3 林芝市公路沿線地質災害易發性區劃

4 地質災害氣象預報模型研究

4.1 地質災害與降水誘發災害概率關系研究

本研究基于前期整理的2010—2019年地質災害資料和降水資料，將地質災害是否發生作為因變量，日降水量作為自變量，以邏輯回歸方法作為耦合手段分析地質災害與降水誘發災害概率關系。林芝市截止至2019年約有66個氣象站點的日降水數據，由于2010—2019年增加或棄用的氣象站點較多，則利用ArcGIS鄰域分析工具，分別提取2010—2019年距歷史災害點最近雨量站的日降水量平均值作為災害點的降水量。相比于降水量空間插值的方式，此方法保證了災害發生當日及前期降水量的準確性。樣本點選擇方式是選擇歷史地質災害點中有準確降水量的樣本點152個，并加入同等數量的無災害發生樣本點，其中，無災害發生樣本點的日期選擇方式是選取2010—2019年汛期(5—9月)較長時間無災害發生的日期，并隨機分配給無災害發生樣本點。

基于災害點樣本、無災害發生樣本和前6日的日降水數據，構建邏輯回歸模型。擬合發現，前4日、前5日和前6日降水量未通過顯著性檢驗。剔除后重新建模，最優結果的Nagelkerke方差可達到0.516，Cox & Snell R方差可達到0.387，各建模因子均通過顯著性檢驗(表3)。

降水誘發地質災害概率的公式為：

(4)

式中：Y為降水誘發地質災害概率值；D0，D1，D2，D3分別為災害發生當日、前1日、前2日、前3日降水量。

表3 研究區建模因子邏輯回歸方程

4.2 地質災害氣象預報模型建立

根據劉傳正等人對國內外地質災害氣象預報模型原理的總結[37]，本研究選取顯式統計預報模型的方法，基于邏輯回歸方法，對地質災害易發性概率H值與降水誘發地質災害概率Y值進行分析，最終得到研究區內地質災害發生的概率值[38]，其計算公式為：

T=a·H+b·Y

(5)

式中：T為地質災害發生概率值；H為地質災害易發性概率值，是地質要素對地質災害發生的綜合體現；Y為降水誘發地質災害概率值；a為易發性區劃占地質災害發生概率權重系數；b為降水因素占地質災害發生概率權重系數。

經擬合，T模型邏輯回歸結果的Nagelkerke方差可達到0.681，Cox & SnellR方差可達到0.51，模型擬合度較高，通過了顯著性檢驗，獲得最終的系數詳見表4。由于模型最終結果是概率值，因此，本研究將表4中的系數歸一化至0～1之內，系數四舍五入后得到a=0.373，b=0.627。

表4 研究區T模型邏輯回歸方程

4.3 地質災害氣象預報模型檢驗

為了確保模型的可操作性，本研究選取部分歷史地質災害點對模型進行精度驗證，將日降水量數據和地質災害易發性概率H值放入地質災害氣象預報T模型中，計算出地質災害發生概率T值。再結合統計分析及專家經驗，確定T值超過0.4發布黃色預警，超過0.6發布橙色預警，超過0.8發布紅色預警。

4.3.1 模型預測準確率計算 選取49個未參與建模的災害點數據作為模型的檢驗數據。將49個災害點對應的當日、前1日、前2日和前3日降水量數據和地質災害易發性概率H值放入地質災害氣象預報T模型中計算，得到不同預報等級的災害點個數，預報結果詳見表5。從表5來看，預報區85.71%的災害點達到黃色預報級別，14.29%的災害點達到橙色預報級別，14.29%的災害點達到紅色預報級別。所以，該模型的預報準確率為85.71%，模型的漏報率為14.29%。

表5 研究區地質災害模型預報效果等級分布

4.3.2 模型個例檢驗 本研究以2020年7月11日318國道林芝縣至波密縣路段多處發生的地質災害為例，驗證地質災害氣象預報模型的應用效果。7月8—11日連續4日318國道林芝縣至波密縣路段有強降水，7月8—11日降水量最大值分別為23.0，35.6，77.3和78.5 mm，將這4日降水量數據帶入地質災害氣象預報模型計算得到2020年7月11日地質災害預報結果(圖4)。據統計，2020年7月11日318國道林芝縣至波密縣路段發生的地質災害點有3個，均達到了橙色預報級別，其中，達到了紅色預報級別的災害點有2個。至此，本研究認為該模型的預報準確性較高，可以為林芝市境內公路沿線地質災害氣象預報提供一定的依據。

圖4 研究區2020年7月11日地質災害預報結果

5 結 論

(1) 針對林芝市公路沿線的地質特性，選取了斷層密度、土地利用、土壤類型、坡度、溝谷密度、年降水量、裸巖率和距河流的距離8個地質災害易發性區劃因子，并通過信息量及邏輯回歸方法確定各區劃因子的系數，系數分別為0.742，0.407，0.871，0.275，0.251，0.427，0.279，0.609。根據林芝市公路地質災害易發性區劃結果圖可以看出，地質災害易發性概率超過0.8的主體區域主要分布在墨脫縣境內扎墨公路、察隅縣東部的察佐公路及察丙路附近；大部分區域易發性概率小于0.4，總體趨勢是離道路越近，易發性概率越大。

(2) 基于地質災害易發性概率值和災害發生當日、前1日、前2日和前3日的日降水數據，通過邏輯回歸方法構建了林芝市公路地質災害氣象預報模型。通過對模型進行檢驗發現，模型預報準確率為85.71%，漏報率為14.29%。

(3) 受林芝市公路沿線地形地貌影響，氣象站點分布不均勻，高海拔地區站點少，地質災害發生地的實際降水與降水資料有偏差。同時，歷史地質災害點信息監測及采集困難，存在錯報及漏報情況。所以，本研究的預報結果可能與實際情況有一定的偏差，但能反映基本趨勢。在未來的研究中，研究人員可通過災害數據篩選、氣象預報模型優化等方式降低模型漏報率，保證實際運行時的預報準確度。