張 聰,姚令侃,2,3,黃藝丹,2,3,邱燕玲,譚 禮
(1.西南交通大學土木工程學院,四川 成都 610031;2.西南交通大學高速鐵路線路工程教育部重點實驗室,四川 成都 610031;3.西南交通大學陸地交通地質災害防治技術國家工程實驗室,四川 成都 610031;4.中鐵二院工程集團有限責任公司,四川 成都 610031)
青藏高原南緣是世界上中低緯度地區最大的現代冰川分布區,發育了世界上罕有的海洋性山谷冰川,并由于冰川的活動形成了大量冰磧湖[1].海洋性冰川具有消融量大、活動性強的特點,使得青藏高原南緣,尤其是喜馬拉雅山脈,成為全球冰磧湖潰決災害事件最密集的地區之一[1-2].冰磧湖位于高海拔地區,一旦潰決,其潰決洪水在重力作用下往往給下游的生命財產和交通設施帶來難以估量的損失.例如,1981年喜馬拉雅中部波曲河流域的次仁瑪錯冰磧湖潰決事件,使得中尼公路超過55 km路段損毀,交通中斷,并導致數百人死亡[1].因此,預測、評估冰磧湖潰決風險對冰磧湖風險調控具有重大的現實和指導意義.
長期以來,隨著國內外學者對冰磧湖潰決災害研究的不斷深入,從潰決的激發因素與潰壩機制出發,對冰磧湖危險性進行評估,取得了一系列成果.徐道明等[3]通過對我國西藏地區潰決冰磧湖的考察,提出了從壩體特征、氣候因素判別危險冰磧湖的幾類指標.呂儒仁等[4]分析了我國西藏冰磧湖潰決特點,統計了有利于潰決事件發生的7類指標.Huggel等[5]利用遙感影像從冰磧湖、壩、激發條件、下游洪水路徑等方面提出冰磧湖潛在危險性評判的三級準則和相應指標.基于上述指標體系可對冰磧湖潰決危險性作出定性評價.
從冰磧湖潰決危險性定量評估角度,呂儒仁等[4]將冰磧湖與母冰川變坡點(或裂隙面)下游冰體的體積比定義為危險性指數,提出了評估由冰崩涌浪造成壩體潰決的定量化指標.陳曉清等[6]從喜馬拉雅地區冰磧湖相關因子中選取7個參數作為直接判別的依據,結合冰磧湖潰決危險指數建立起一種定性-半定量的個體冰磧湖危險性評價方法.舒有鋒等[7]從冰川、湖盆、終磧堤、氣候等特征中提取冰磧湖潰決危險性評價指標,通過粗糙集理論確定各指標權重,建立了冰磧湖潰決危險性評價方法.王欣等[8]將冰磧湖潰決機制的定性描述與潰決概率進行轉換,建立了冰磧湖潰決的事故樹模型,對我國青藏高原以南冰磧湖潰決危險性進行了分析.
冰磧湖潰決通常是一種因素主導,多種因素共同作用的具有高度不確定性的復雜物理過程,不同誘因組合可能導致不同潰壩模式以及嚴重程度不同的潰決洪水災害.上述評價方法對冰磧湖潰決這類高度不確定性、多狀態、共因導致的失效問題針對性不強.貝葉斯網絡(B-Ns)分析法具有良好的適用性[9-10],能夠很好地彌補以上評價方法的不足,同時可充分利用已有案例作為先驗知識,提高預測結果的客觀性和準確性,但目前未見其應用于冰磧湖潰決危險性分析.
中巴鐵路、中尼鐵路等穿越冰磧湖分布區的口岸鐵路均列入國際通道的規劃中.這些跨境鐵路將穿越喜馬拉雅山脈,冰磧湖潰決已成為對鐵路線路方案起控制作用的重大災害類型,而冰川覆蓋區是我國鐵路沒有建設經驗的地貌單元.冰磧湖一般位于極高海拔區,難以通過工程手段處治,因此以在選線階段開展風險調控最為主動.鑒此,以冰磧湖潰決危險區鐵路方案風險評價為目標,首先對青藏高原南緣74例冰磧湖潰決災害實例潰決機制進行統計分析;在此基礎上針對冰磧湖潰決影響因素眾多、關系復雜且具有高度不確定性的特點,利用多態BNs建立了冰磧湖潰決概率預測模型;依據冰磧湖潰決危險性評價結果與沿河線河谷地貌形態特征提出了針對冰磧湖潰決危險性的鐵路選線方案評價方法;最后,以中尼鐵路跨喜馬拉雅段吉隆、樟木局部走向方案比選為例說明評價方法的作業程式.為冰磧湖潰決危險區鐵路方案風險評價提供具普適性的技術評價體系.
B-Ns能夠有效地在不確定性環境中實現概率推算與結果預測,利用其進行危險性分析的一般流程是構建B-Ns模型、確定節點變量及相應的概率分布,然后進行推理計算.其中節點變量及網絡構建均依賴于已有案例作為先驗知識,因此,本節首先對青藏高原南緣冰磧湖潰決災害案例進行梳理.
結合實地考察、影像資料與歷史資料[11-17],本文共梳理了青藏高原南緣74例冰磧湖潰決災害事件(圖1,附加材料表S1).可知,該區域冰磧湖潰決事件主要發生在喜馬拉雅中、東部以及藏東南地區.表S1中,對冰磧湖潰決機制有記載的案例為62例,統計知其潰決機制可以大致分為5類(表1),其中冰崩涌浪漫頂是最主要潰決機制.此外,有8例為多種潰決機制共同作用導致.

圖1 青藏高原南緣冰川、冰湖及冰磧湖潰決災害分布Fig.1 Distributions of glaciers, glacier lakes and morainedammed lake outburst events in the southern Tibet Plateau
1.2.1 構建多態B-Ns模型
構建B-Ns模型的首要工作是確定B-Ns節點變量.將表1中冰磧壩失效的5種主要模式作為一級指標,然后結合冰磧湖潰決災害形成的基本條件,分別對5個一級指標建立各自的二級分析指標.

表1 青藏高原南緣冰磧湖主要潰決機制統計Tab.1 Main outburst mechanisms of the moraine-dammed lakes in the southern Tibet Plateau
1)冰崩涌浪漫頂(y1):冰川末端的冰體在氣候等激發下,突然脫離母冰川進入下游冰磧湖激發冰磧湖涌浪,同時導致冰磧湖水位驟漲,漫頂沖刷冰磧壩致使壩體潰決.冰崩的斷裂面一般位于冰川坡度突變處或者冰舌裂隙發育部位,因此,該部位下游冰體可稱為危險冰體.y1的影響因子分析如下:
① 危險冰體坡度(x1)決定著冰體崩落起動的難易程度.由青藏高原南緣的冰磧湖潰決案例統計可知,發生冰崩災害的冰川末端坡度變化較大,但當冰川末端的危險冰體坡度超過8° 時,容易導致冰崩災害的發生[6].
② 危險冰體指數(x2)為危險冰體體積與冰磧湖湖水體積的比值.危險冰體滑入冰磧湖后產生涌浪,同時導致冰磧湖水位驟漲.二者綜合形成的水頭高度越大,壩體潰口處起動泥沙顆粒也越容易起動,相應發生潰決可能性越高[18].因此可以采用危險冰體指數表示其危險性.參照冰磧湖潰決案例,潰決前危險冰體指數一般位于0.11 ~ 0.52[6].
③ 冰舌冰磧湖距離(x3).冰崩涌浪的首要形成條件是冰川冰磧湖磧處于鄰近狀態.參照冰磧湖潰決案例,已潰決的冰磧湖x3值大多小于1 km.
④ 冰川裂隙發育程度(x4).在氣候反常年份,冰川容易在裂隙發育面發生斷裂,因此冰川裂隙發育程度與冰崩災害的易發性呈現正相關關系.
⑤ 年際氣溫變異系數(x5).氣溫突增使得冰川底部融水增加,摩阻減小,同時導致冰川內部應力差增加,容易引發冰崩.采用年際氣溫變異系數代表氣溫突變程度,系數值越大,氣溫年際變化越劇烈,也越有利于災害的形成.由西藏地區冰磧湖潰決區鄰近氣候站溫度統計知,冰磧湖潰決區的氣溫變異系數值位于0.23 ~ 0.27.
2)管涌(y2):壩體厚度和壩體物質組成是管涌的影響因素[1].目前,對已潰決冰磧湖統計發現大部分冰磧壩頂寬(x6)不超過60 m.其次,堰塞壩體的物質組成影響著滲流穩定,特別是當冰磧壩體內部含有埋藏冰時,隨著氣溫升高,埋藏冰融化容易引發冰磧壩的管涌現象.當前并無判別冰磧壩內是否存在埋藏冰的方法,但是徐道明等[3]調查后認為,在小冰期末時期所形成的冰磧堰塞湖的壩體中,埋藏冰是普遍存在的,且受氣溫變化的影響較為明顯.據此,可以采用冰磧壩頂寬、壩體物質狀態(x7)與年際氣溫變異系數對管涌的危險性進行分析.
3)洪水漫頂(y3):在氣候異常年份,冰川融水與降雨導致的超標洪水對冰磧湖補給超過溢流、蒸發的損耗時,冰磧湖容易發生漫流型潰決.采用年際降雨變異系數表征降雨的不均勻度,則冰磧湖漫流潰決影響因子可以采用冰磧湖滿溢程度(x8)、年際降雨變差系數(x9)與年際氣溫變異系數表示.
4)冰磧壩坍塌(y4):冰磧壩多為松散巖屑或積雪(冰)夾帶巖屑,穩定性較差的自然堆積體.壩體的坡度越大,其穩定性越差,在內、外部因素作用下容易坍塌[1].對青藏高原南緣冰磧湖潰決案例調查表明,潰決冰磧堰塞壩背水面的坡度多位于20° ~ 33°.其次,壩頂若存在漫流下切槽道,則壩體容易在溯源侵蝕的作用下坍塌.此外,冰磧壩內埋藏冰消融也可能導致壩體坍塌.因此,可以采用背水坡度(x10)與壩頂漫流槽道分布狀態(x11)以及年際氣溫變異系數分析.
5)外界動力激發(y5):外界動力主要為地震作用.地震力一方面促使冰磧湖附近巖體、冰雪崩滑激發涌浪,另一方面地震直接激發水體涌浪導致漫頂現象的發生.Ai等[18]利用地震振動臺模擬地震與滑坡共同作用下復合涌浪實驗發現,地震峰值加速度(PGA,x12)與涌浪波高正相關,當x12> 0.25g(超越概率等于10%)時,涌浪波高可達到冰磧壩顆粒起動臨界條件,可能引發冰磧湖潰決災害.
確定上述5種潰壩機制作為中間節點變量,根據節點變量之間的因果關系建立B-Ns拓撲結構,如圖2所示,其中:T為葉節點事件,即冰磧堰塞湖潰決;yj為中間節點事件 (j=1,2,···,m);xi為根節點事件 (i=1,2,···,n).依據有利于形成災害的節點狀態分界值,對B-Ns根節點變量進行離散化處理,得到各根節點變量狀態及取值范圍,如表2所示.表中:狀態值0、1分別表示安全、危險狀態.

圖2 冰磧湖潰決B-Ns模型Fig.2 B-Ns model of moraine-dammedlakes outburstrisk

表2 根節點因素多區間狀態劃分標準Tab.2 Division of risk factors in multi-interval safety states for root nodes
傳統的事故樹模型以及二態B-Ns系統僅將節點狀態簡化為“正常”和“失效”兩種,但在實際中中間節點與葉節點往往存在著中間狀態,即其節點狀態雖然不屬于正常狀態范疇,但是并未導致大的災害形成.以冰磧湖潰壩機制中占比最高的冰崩涌浪漫頂為例,漫頂涌浪往往需要達到一定的水頭高度才能使得堰塞壩潰口處粗顆粒泥沙連續起動形成潰決連鎖反應[18].因此,根據事故嚴重程度,本文將中間節點y1的狀態分為低、中、高危險3類,狀態值分別用0、1、2表示;其余潰決機制(y2~y5)由于發生幾率較小,為降低貝葉斯網絡復雜度,仍將其簡化為兩個狀態(以0、1表示).對于葉節點T,從堰塞湖潰決物理過程角度,其潰決方式可以分為全潰與部分潰,其中全潰形成的相對洪峰流量最大,危險性也最高[19].根據潰決危險性嚴重程度,將葉節點T分為低、中、高3種危險狀態.
1.2.2 節點概率分布的獲取
構建B-Ns結構后,需要確定根節點各狀態先驗概率分布及其他節點(中間節點和葉節點)條件概率分布.節點概率的獲取方法主要有資料統計法和專家知識法兩種.B-Ns方法優勢之一在于能夠有效地融合多源信息,既可以利用現有資料統計獲取客觀概率值,也能夠利用專家知識對節點概率進行補充.因此,本文在確定B-Ns節點概率分布時結合資料統計法與專家知識法互相驗證.
根節點先驗概率可依據資料統計法獲得.根據所研究區域的冰川、冰磧湖、遙感影像與數字高程模型(DEM)地形對根節點因子落入各狀態區間頻數進行統計,從而確定根節點各狀態的先驗概率.本文所使用統計資料中,冰川、冰湖數據來自國際山地綜合發展中心(ICIMOD,2015年)與中科院寒旱所(2014年)提供的第二次冰川編目;冰磧壩等幾何與結構參數源自對遙感影像(2014年—2019年)與DEM地形(2009年)的解譯;氣候數據來自中國氣象數據網提供的研究區氣象站近30年氣象資料.冰磧湖潰決案例為本文的附加材料表S1.
通過以上方法獲取根節點先驗概率之后,仍然需要確定其他節點的條件概率表(CPT).節點處CPT可利用模糊數學法將專家的定性描述轉化為定量的發生概率[8].針對專家知識法獲得的CPT存在主觀性的問題,可以將獲取的條件概率與實際案例互相驗證,以降低主觀性,操作程式見1.3節案例.
1.2.3 冰磧湖潰決概率計算
根據B-Ns的運算規則,可以由根節點的先驗概率分布出發正向推算葉節點事件的發生概率.
在多態B-Ns中,根節點為xi;中間節點為yj;葉節點為T.、和Tq分別代表根節點、中間節點、葉節點的風險狀態.其中:ai=0,1,···,ui?1;bj=0,1,···,vj?1;q=0,1,···,r?1,ui、vj、r分別為根節點、中間節點、葉節點的狀態數.則當根節點各狀態的概率分別為P,P,···,P時,葉節點T處于風險狀態Tq的概率計算式為

式中:π(yj)為中間節點yj的父節點集合;π(T)為葉節點T的父節點集合.
依據待評估冰磧湖的證據信息,利用B-Ns的正向推理可以計算葉節點T各狀態的發生概率.確定各風險狀態的隸屬度函數為

式中:p為葉節點各狀態發生概率值.
根據最大隸屬度原則確定葉節點隸屬的狀態.評估程式如圖3所示.

圖3 基于B-Ns的冰湖潰決危險性評估程式Fig.3 Risk assessment procedure for moraine-dammed lake outburst based on B-Ns method
評價模型的精度直接關系到評價結果的可靠性,因此使用前有必要對該模型進行檢驗.本節利用1981年7月11日波曲河流域次仁瑪錯冰磧湖案例對B-Ns模型的精度與作業流程進行驗證和說明.
根據呂儒仁等[4]的調查,潰決當年喜馬拉雅南坡的氣候背景為暖干氣候,降水因為量小而作用很小,突然升溫導致冰崩,疊加管涌效應導致冰磧湖潰決.該冰磧湖的海拔高度為4660 m,潰決前夕冰舌到達海拔4700 m位置,距離冰磧湖非常近.冰川前緣地形為陡崖狀,且在陡崖處呈不連續狀態,裂紋發育.潰決之前,冰磧湖處于滿溢狀態,壩頂槽道發育,已有湖水從壩體中連續滲出.潰決危險性評價流程如下:
1)確定節點概率分布
根節點各風險狀態區間的先驗概率可依據資料統計法獲得,如表3所示.

表3 冰磧湖潰決根節點的先驗概率Tab.3 Prior probabilities of root nodes of moraine-dammed lake outburst
依據各風險因素間的因果關系,結合案例資料與專家知識可構造中間節點y1~y5及葉節點T各狀態區間的條件概率分布.本文僅選取中間節點y1部分條件概率分布作為示例,如表4所示.表中,第一行數字的意義為:當根節點x1~x5狀態量均為0時,中間節點y1狀態量為0的概率等于1.0,為1、2的概率均為0.其余各行數字含義可類似推出.完整表格詳見附加材料表S2.
將表3中根節點x1~x5的先驗概率與表4中y1的條件概率代入式(1),可求得節點y1各風險狀態的概率為同理可得,P(y1=1)=0.2655,P(y1=2)=0.0526.

表4 中間節點y1的條件概率Tab.4 Conditional probabilities of node y1

上述計算結果表明,在沒有相關證據信息輸入B-Ns模型條件下,中間節點y1的狀態為低危險的概率很大,中、高危險概率均較小,其中為高危險的概率僅為0.0526.據統計,波曲河流域冰湖數量為121個,參照附加材料表S1,該流域發生嚴重冰崩涌浪漫頂潰決災害的案例為4次,概率值計算為0.0331,與本文計算的P(y1=2)=0.0526 較為接近.可初步驗證B-Ns模型結構及參數.根據類似方法,可構造其余節點的條件概率.
2)依據證據信息對待評估冰磧湖潰決危險性進行評價.
根據潰決前次仁瑪錯冰磧湖狀態,更新貝葉斯網絡中各根節點狀態量,結合中間節點的條件概率進行推理計算,可得P(T=0)=0.0670,P(T=1)=0.2381,P(T=2)=0.6949.根據最大隸屬度原則,該冰磧湖的潰決危險性評價結果為高危險,與實際情況相一致,從而對本文B-Ns模型準確性進行了進一步驗證.
各類風險區線路長度應是潰決危險區鐵路選線評價時最重要的指標,其計算需綜合考慮冰磧湖的危險性與線路工程所處河谷段地形.
青藏高原地區,冰磧堰塞湖均分布于主河兩側的支溝內,而線路工程多沿主河布設.冰磧湖一旦潰決,潰決洪水先沿著支溝演進,進入主河后若形成超標洪水才會對線路工程造成災害.能否形成超標洪水的必要條件是線路工程位于潰決洪水的淹沒范圍內,而每個冰磧湖的淹沒范圍與其體積(V)有關.現借用Washakh等[20]提出的冰磧湖潰決洪水淹沒長度(L)計算式 (L=0.52V+36.13),若定義潰決洪水到線路工程的路徑長度為Ld,則篩選出需進行風險評估的冰磧湖篩選判據為

能夠形成超標洪水的另一必要條件是沿河線所在河谷地形,顯然同樣流量的條件下越是地形狹窄的河道形成超標洪水的情況越嚴重.
綜上,各類風險區線路方案長度指標計算不僅需要考慮冰磧湖潰決危險性,還需要考慮冰磧湖潰決影響范圍、沿河線河谷形態等因素,具體計算流程如圖4所示.

圖4 各類風險區線路長度計算流程Fig.4 Calculation process of the railway line length in different risk areas
鐵路方案冰磧湖潰決危險性評價不僅需要考慮流域內冰磧湖的自然狀態(數量、總面積、歷史潰決次數),也要考慮危險性冰磧湖對鐵路線路的影響.結合這兩類考慮因素,提出鐵路方案冰磧湖潰決危險性評價指標體系,如表5所示.以中尼鐵路為例進行作業程式說明.

表5 鐵路方案冰磧湖潰決危險性評價指標Tab.5 Evaluation indexes of railway schemes considering moraine-dammed lake outburst risk
中尼鐵路是擬建穿越喜馬拉雅山脈的鐵路,線路從日喀則出發先在高原面上行走,而后下高原面至加德滿都,這時面臨經吉隆廊道與經樟木廊道這兩個局部走向方案的比選問題(圖5).吉隆、樟木廊道均位于喜馬拉雅中部希夏邦馬峰冰川群分布區,冰磧湖數量眾多,冰磧湖潰決洪水已成為對線路方案起控制作用的重大災害類型.

圖5 吉隆、波曲流域危險冰磧湖與線路方案示意Fig.5 Railway schemes and distribution of dangerous moraine-dammed lakes in the Gyirong and Poiqu river
1)研究區自然條件
吉隆流域位于喜馬拉雅中部.東西走向的喜馬拉雅山脈將吉隆流域分割成南北兩部分,北部海拔較高、地勢平緩;南部海拔低,但是山體陡峻,高差較大.遠古時期吉隆地區分布有吉隆、沃瑪等湖泊,由于青藏高原的不均勻抬升,導致吉隆藏布溯源侵蝕切穿原有湖泊,湖水下泄,這樣原古湖盆與現代河流形成了吉隆藏布上游段寬窄相間的河谷地貌形態.此外,由于海洋性冰川塑造地貌效應,中游段也形成了一些寬谷地形.寬谷和峽谷段分布如圖5所示.
樟木廊道位于波曲流域,波曲河發源于西藏聶拉木縣波絨鄉,經聶拉木縣城到樟木鎮.流域內地形高差較大,最高峰為海拔超過8 km的希夏邦馬峰.波曲的河谷地貌以聶拉木縣城為分界線,其北基本位于高原面上,主要以寬谷地貌為主,其南屬于溯源侵蝕河段,以峽谷地貌為主(圖5).
2)研究區冰磧湖潰決危險性評價
首先利用式(2)篩選出需要進行危險性評價的冰磧湖,然后利用本文建立的B-Ns評價模型,對吉隆、樟木廊道流域內冰磧湖潰決危險性進行評價(如圖5).可知:吉隆流域危險性較大的冰湖總數為7個(其中高危險冰磧湖4個,中危險冰磧湖3個);樟木流域危險性較大的冰磧湖總數為13個(其中高危險冰磧湖6個,中危險冰磧湖7個).
3)中尼鐵路樟木、吉隆局部走向方案比選
中尼鐵路跨喜馬拉雅段的樟木、吉隆局部走向方案如圖5所示.
依據本文的鐵路方案冰磧湖潰決危險性評價指標對兩方案進行比選,各項指標值如表6所示.由表可知:雖然樟木方案在冰磧湖分布區的線路長度小于吉隆方案,從選線角度屬于線路短直方案,但是樟木方案的高、次高、中風險區段線路長度都遠超過吉隆方案.因此,從冰磧湖潰決風險角度,吉隆方案優于樟木方案.

表6 吉隆、樟木線路方案比選表Tab.6 Comparison of the railway schemes in the Gyirong and Poiqu river basins
本文以冰磧湖潰決危險區鐵路方案風險評價為目標,建立了基于多態B-Ns的冰磧湖潰決危險性分析方法及冰磧湖潰決風險區新建鐵路選線方案風險評價指標體系.主要結論有:
1)通過對青藏高原南緣74個冰磧湖潰決災害實例統計分析,在此基礎上針對冰磧湖潰決影響因素眾多、關系復雜且具有高度不確定性的特點,建立了多態B-Ns冰磧湖潰決災害發生的概率預測模型;以冰磧湖分布區沿河線路工程為承災體,依據冰磧湖潰決危險性評價結果與沿河線河谷地貌形態特征,提出了鐵路選線方案冰磧湖潰決危險性評價指標體系,從而為冰磧湖潰決危險區鐵路方案風險評價提供了具普適性的技術評價體系.
2)本文提出的鐵路方案冰磧湖潰決危險性評價指標,僅是從冰磧湖角度對現行評價體系進行增補.在實際工程中,仍需在確定常規方案技術經濟指標的基礎上,結合冰磧湖潰決危險性評價指標進行綜合分析,評選最合理方案.
3)鐵路選線一般需經過從線路走向的擬定到空間定線等從粗到細逐步優化的階段,直至確定合理的線路位置.本文提出的鐵路選線方案冰磧湖潰決危險性評價指標體系主要用于鐵路行經地區局部走向方案的選擇;此外針對空間定線作業階段,提出在冰磧湖潰決高危險區的峽谷段,宜采取抬高程、盡量不跨河、酌情預留采取提高限坡措施條件等定線要點;二者共同構成冰磧堰塞湖潰決減災選線風險調控技術體系.
備注:附加材料在中國知網本文的詳情頁中獲取.
致謝:本文研究得到了中國科學院水利部成都山地災害與環境研究所蘇立君研究員,中鐵第一勘察設計院集團有限公司苗曉岐教授級高工,梁棟、孫先鋒高級工程師等熱情的幫助和指教,在此一并深致謝意;感謝中國科學院國際伙伴計劃“一帶一路”科技合作項目“中尼交通廊道災害風險及其應對策略研究”(131551KYSB20180042)、中鐵二院工程集團有限責任公司科研項目(KYY2020054(20-22)對本文研究的資助.