Primary establishment of an early warning model of debris flow hazards in Nyingchi City of Tibetan autonomous region based on raster runoff simulation
-
摘要: 西藏林芝市泥石流灾害频发，亟需建立泥石流灾害预警模型，预测林芝市泥石流灾害可能发生的区域，减少泥石流灾害导致的损失。文章提出了一种基于栅格径流汇流的林芝市泥石流灾害预警模型，从栅格像元尺度上模拟流域各位置上的水深，以提高泥石流预警的空间针对性。该模型将泥石流致灾因子分为背景因子和激发因子。通过林芝市裸岩率、河床纵比降等因子的逻辑回归，获取林芝市泥石流灾害概率，作为泥石流预警模型的背景因子；引入栅格径流汇流模型，以站点降水和雪水当量为模型的水量输入，模拟预警时段内的流域各位置上的模型水深，作为泥石流预警模型的激发因子。利用二元逻辑回归的方法计算背景因子和激发因子的权重，建立泥石流预警模型。利用2011—2020年18次历史灾害对模型进行验证，落入预警区内的灾害点占比64.4%，预警精度较高，对于林芝市泥石流灾害预警具有一定的指导意义。Abstract: Debris flow disasters occur frequently in Nyingchi City, Tibet. There is an urgent need to establish an early warning model to predict the possible areas of debris flow disasters in Nyingchi City and reduce the losses caused by those disasters. This paper presented an early warning model based on raster runoff simulation in Nyingchi City, which can simulate the water depth at each location in the watershed and improve the spatial pertinence of debris flow early warning. In this model, the disaster factors of debris flow are divided into background factor and excitation factor. The probability of debris flow disaster in Nyingchi City is obtained by logistic regression of many factors such as bare rock rate, vertical slope of riverbed and so on, which is used as the background factor. The raster runoff simulation model is introduced to simulate the predicted water depth at each position of the basin during the early warning period by importing the precipitation and snow water equivalent data. Using binary logistic regression method to calculate the weight of background factor and excitation factor, the final model is obtained. 18 historical disasters from 2011 to 2020 were used to verify the model. The disaster points falling into the early warning area accounted for 64.4%, which has certain guiding significance for the early warning of debris flow disaster in Nyingchi City.
-
Keywords:
- Nyingchi City /
- debris flow disaster /
- runoff simulation
-
0. 引言
西藏林芝市地势海拔高，多高山峡谷，冰川覆盖度高，是全国泥石流灾害主要易发和高发区，尤其在4—9月份雨季期间，林芝市地区降水量充足，极易在当地地形地质条件下产生泥石流灾害，对经济和人口安全造成危害；另一方面，由于全球变暖，林芝市冰川发生大面积消融，林芝市七个县（巴宜区、米林、工布江达、朗县、波密、察隅；墨脱县无资料）每年都要发生几十次程度不同的泥石流灾害，造成人员、牲畜伤亡以及房屋、公路冲毁等重大损失，泥石流灾害成为严重制约林芝市经济社会发展的问题[1-3]。
目前泥石流灾害预警模型主要分两类，一类是仪器接触型，通过分析泥石流灾害发生机理以及环境数据建立预警模型，二类是非接触型，通过泥石流灾害因子统计分析后进行泥石流灾害预警[4-6]。非接触式分析法主要是通过层次分析法、熵值法、灰色理论分析法、模糊数学分析法等等统计学方法分析与降雨量之间的关系，建立泥石流灾害预警模型[7]，是目前泥石流面域预警的主流方法。王英杰等[8]从岷江上游流域的地质环境背景出发，分析了泥石流分布情况、灾害特征以及所处的发育程度，在划分子流域的基础上，通过采用灰色系统理论和确定降雨临界值建立了岷江上游预警报模型；M Ponziani等[9]根据水文模型和降雨阈值发布水文地质警报，以确定该地区不同地区触发浅层滑坡的可能性；M Martinengo等[10]基于反向动力学方法（BDA）确定石质泥石流降雨阈值的不确定性，用于估算与实测泥石流密切相关的降雨持续时间和平均强度。
由于临界雨量法以流域整体为单元进行预警，无法预测泥石流灾害发生部位。泥石流本质上也是一种流体运动，与流域径流存在较大的联系。近年来，将降雨—径流—土壤渗透相结合，从泥石流产生过程构建泥石流灾害预警模型成为研究热点。包红军等[11]采用GBHM分布式水文模型模拟泥石流发生的水文过程；基于无限边坡模型，利用安全系数判别坡面失稳过程；以流域径流和失稳土体耦合后的水土混合物密度作为泥石流预报标准实现区域泥石流预报；王治华等[12]研究了流域内泥石流的形成和传播过程，将其分为降雨入渗、径流和泥石流演进三个阶段，建立了一个综合考虑初始条件、运动机制和卷吸效应的径流型泥石流模型，有效地描述了泥石流形成和传播过程中各阶段的行为特征。
目前，多数径流汇流采用（半）分布式水文模型，一般以流域为单位，仅能获得流域出口的流量过程[13-15]。本文尝试提出一种栅格模式的径流汇流模拟，从栅格像元尺度上模拟流域各位置上的水深，以模拟水深替代传统的面雨量，并综合沟口距离因子、危险性概率建立泥石流灾害预警预报模型，以提高林芝市泥石流灾害预警的空间针对性，为林芝市防灾减灾提供科学依据。
1. 研究区概况与模型构建思路
1.1 研究区概况
林芝市地处西藏东南部，平均海拔为2700 m左右，位于喜马拉雅山脉附近的峡谷部分属高山峡谷区，山高谷深。林芝市泥石流沟都具有完整而典型的山谷型泥石流流域特征，岭谷高差为1384～4596 m，山体坡度为23°～38°，沟床比降大都在25%～55%，赋予了泥石流灾害有利的地形条件。区内土壤多半淋溶土、半水成土及淋溶土，为泥石流提供黏粒部分。在湿润气候下极易发生黏土化，通过多种形式可为泥石流灾害提供大量物源；同时，林芝市多断层结构岩质疏松，在风化作用或者雨水的冲刷下极易产生松散固体物质并且发生移动，为林芝市泥石流灾害的产生提供了良好了物质基础；另外，该区域较大部分地区是第四纪冰川作用盛行地带，除风化，构造破碎形成的碎屑物质外，分布最广最多的是冰磧物质和冰水沉积物。已有研究资料表明，在古乡沟中上游段内的围谷盆地中就堆积储量约为4.00×108 m3的新老冰碛物，是该沟发生大规模泥石流并频繁暴发的主要因素之一。因此，林芝市地质环境脆弱，地形起伏大，土质基本为沙壤结构疏松易结团，是泥石流灾害频繁发生的重要因素。
从空间分布来看，泥石流灾害主要集中在波密县、察隅县、林芝县、墨脱县。其中波密县、察隅县发生泥石流灾害频率最高是重点灾害防护地区，占历史灾害总频率58%，米林县、朗县、工布江达县发生灾害频次较低（图1）。
![]()
图 1 林芝市泥石流灾害点、隐患点分布图Figure 1. Distribution of debris flow disaster sites and hidden danger points in Nyingchi City1.2 泥石流灾害模型构建思路
泥石流灾害是地质、地形、地貌、土壤、植被、降水与冰川融水等因素共同作用的结果。相等规模的降水或融雪，因为地形、地质土壤等条件不同而产生不同程度的泥石流灾害。由于地质、地形、地貌等因素在一定的时间内保持相对的稳定，将其作为区域的背景因子，降水和融雪作为泥石流灾害的激发因子，构建泥石流灾害预警模型。
首先，对林芝市进行小流域划分。将背景因子纳入到不同的流域中，分析流域的背景因子对泥石流灾害的影响，构建基于流域单元的林芝市泥石流灾害的危险性概率模型；然后，引入降雨融雪—径流汇流模型模拟径流水深，将模拟水深作为泥石流灾害激发因子；最后，建立基于径流模拟水深、危险性概率背景值因子、沟口距离因子的林芝市泥石流灾害预警模型。泥石流灾害预警模型构建思路如图2所示。
2. 泥石流灾害危险性概率计算
2.1 林芝市流域划分
目前，泥石流灾害的危险性区划大多建立在单个栅格单元的基础上。单个栅格单元难以量化流域整体特征[16]，如流域面积、纵比降等，而这些因子是构建泥石流模型的重要因子。因此，本文利用区域DEM，基于盆域模型将林芝市划分为大量的小流域，作为泥石流危险性概率计算的基本单元。由于部分流域较细碎手动修正，共计16648个小流域（图3）。
2.2 泥石流灾害诱发因子的选取与逻辑回归分析
诱发泥石流灾害的因素较多，形成机理较为复杂。依据泥石流发生机理，参考前人经验，选择裸岩率[17]、流域面积、纵比降、河流、道路、断层密度、土地利用、土壤类型、隐患点密度、沟谷密度、年降水量因子作为泥石流灾害诱发因子，基于信息量结合逻辑回归的方法确定泥石流灾害危险性概率计算。
各单因子信息量
$I(L, x_{n})$ 的计算方法如下[18-19]：$$ I(L,{x_{{n}}}) = \ln \frac{{{N_i}/N}}{{{S _i}/S}} $$ (1) 式中：
$ S $ ——林芝市评价单元总数；$ N $ ——林芝市地区发生泥石流灾害所分布的总单元数；$ {S_i} $ ——林芝市内含有评价因子$ {x_n} $ 的单元数；$ {N_i} $ ——分布在因子$ {x_n} $ 内特定类别的泥石流灾害单元数。$$ Z = \ln \left(\frac{p}{{1 - p}}\right) = A + {B_1}{X_1} + {B_2}{X_2} + \cdots + {B_n}{X_n} $$ (2) $$ P = \frac{{{\rm{exp}}(Z)}}{{1 + {\rm{exp}}(Z)}} $$ (3) 式中：
$ p $ ——泥石流灾害危险性概率，取值在[0，1]之间；(
$ 1-p $ )——泥石流灾害不发生的概率；$X_{1} ,\cdots. X_{n}$ ——影响地质灾害发生的因子；$B_{1},\cdots ,B_{n}$ ——各影响因子对应的逻辑回归系数。将剔除误差数据后的701个隐患点数据记录作为建模数据，128个有效灾害点数据作为检验数据，提取地质隐患点和未发生地质灾害点致灾因子信息量值，以泥石流灾害是否发生为因变量，裸岩率、流域面积、纵比降、河流、道路、断层密度、土地利用、土壤类型、隐患点密度、沟谷密度、年降水量信息量值为自变量，通过逻辑回归计算得到表1。
表 1 逻辑回归分析结果Table 1. Results of logistic regression analysis指标因子 B S.E Wals df Sig. exp（B） 裸岩信息 0.925 0.456 4.119 1 0.042 2.523 流域面积 0.880 0.394 5.000 1 0.025 2.412 沟床纵比降 1.042 0.186 31.348 1 0 2.834 河流 1.000 0.122 66.763 1 0 2.717 道路 0.879 0.105 69.489 1 0 2.409 断层密度 0.995 0.442 5.059 1 0.024 2.704 土地利用 0.541 0.245 4.888 1 0.027 1.718 土壤类型 0.603 0.112 18.239 1 0 1.828 隐患点密度 0.896 0.297 63.593 1 0 2.449 沟谷密度 1.295 0.211 19.043 1 0 3.650 年降水量 1.064 0.281 25.459 1 0 2.897 常量 −1.485 0.281 27.852 1 0 0.226 注：B为逻辑回归系数；S.E.为标准误差；Wals为卡方值统计量；df为自由度；Sig.为显著性。 各因子均通过显著水平为0.05的Wald检验，且因子间相关系数很小，相互独立；多次随机生成未发生泥石流灾害点并循环建模，模型结果稳定；最终确定泥石流灾害危险性概率如下：
$$ P = \frac{\begin{split} &\exp({{0}}{{.925X_1}} + {{0}}{{.880X_2}} + {{1}}.042{{X_3}} + {{1X_4}} + {{0}}{{.879X_5}} + 0.995{{X_6}} + \\ &{{0}}{{.541}}{{{X}}_{{7}}} + {{0}}.603{{X_8}} + {{0}}{{.896X_9}} + {{1}}.295{{}}{{{X}}_{10}} + 1.064{{{X}}_{11}}{{ - 1}}{{.485}}) \end{split}} {\begin{split} &1 + \exp({{0}}{{.925X_1}} + {{0}}{{.880X_2}} + {{1}}.042{{X_3}} + {{1X_4}} + {{0}}{{.879X_5}} + 0.995{{X_6}} + \\ &{{0}}{{.541}}{{{X_7}}} + {{0}}.603{{X_8}} + {{0}}{{.896X_9}} + {{1}}.295{{X_10}} + 1.064{{X_11}}{{ - 1}}{{.485}}) \end{split}} $$ (4) 式4中，
$ P $ 为泥石流灾害危险性概率值，$ {X}_{1} $ ～$ {X}_{11} $ 分别为裸岩率、流域面积、纵比降、河流、道路、断层密度、土地利用、土壤类型、隐患点密度、沟谷密度、年降水量的信息量值。2.3 泥石流灾害危险性概率背景值的栅格化
由于泥石流灾害危险性概率以流域为单元，仅能得到各流域的危险性值。在本文的预警模型中，径流汇流模拟和泥石流预警预报均是以栅格像元为单元进行的。因此，需要按照统一的空间分辨率，将流域的危险性概率值分配到林芝市250 m的栅格单元中。具体的做法是：首先依据式（4）计算每一个小流域的危险性概率值，然后在250 m分辨率的栅格体系中，将该流域的所有栅格设置为流域的危险性概率值。
2.4 泥石流灾害危险性概率精度验证
将泥石流灾害概率分为五级，利用128个林芝市历史泥石流灾害点来对泥石流灾害概率进行精度验证，各灾害点所处泥石流灾害危险性概率所占比例如表2所示。
表 2 泥石流灾害的概率占比Table 2. Probability proportion of debris flow disaster级别 取值区间 灾害点数量 百分比/% 一级 [0, 0.2) 8 6.25 二级 [0.2, 0.4) 5 3.90 三级 [0.4, 0.6) 2 1.56 四级 [0.6, 0.8) 9 7.03 五级 [0.8, 1) 104 81.25 在五级内的历史灾害点共104个，占总历史灾害点81.25%；四级内的灾害点共9个，占7.03%；三级内的灾害点共2个，占1.56%；在二级及一级内的灾害点共13个，占比10.15%。从表中可见，验证数据与林芝市泥石流灾害危险性概率结果匹配程度较高，为后续泥石流预警模型的建立奠定了基础。
3. 栅格模式下的径流汇流模拟
坡面流和河道径流是泥石流灾害产生的重要条件。提出一种栅格模式下的径流汇流模拟方法，将坡面流和河道径流统一起来，按平流项、压力项和外力项分别进行时间片上每一个栅格的水深和水速计算。通过时间片的迭代，模拟流域的径流汇流过程。
3.1 栅格的水量平衡公式
对于流域的每一个栅格c，建立水量平衡公式：
$$ W_c^{n + 1} = W_c^n + R_c^n + O_c^n + In_c^n + I_c^n + E_c^n + Out_c^n $$ (5) 式中：
$ W_c^{n + 1} $ 、$ W_c^n $ —时间片上中心栅格c输出水深和初 始水深；$ R_c^n $ 、$ O_c^n $ 、$ I_c^n $ 、$ E_c^n $ ——中心栅格时间片上的降水量、冰 川消融量、下渗量、蒸发量；$ In_c^n $ 、$ Out_c^n $ ——模型中平流项；n——时间片起点序号；
n+1——时间片终点序号。
其中，降水量由气象站观测数据插值得到栅格雨量；冰川消融量通过度日模型计算；由于泥石流预警的时段较短，将蒸发率设定为常数；水体在地表运动时，下渗程度受多种因素影响，利用利普下渗曲线进行模拟计算：
$$ f_t^{} = f\left( t \right) = 0.5 \times S \cdot {t^{ - 0.5}} + A $$ (6) 式中：
$ S $ 、$ A $ ——与土壤性质有关的参数。$ In_c^n $ 、$ Out_c^n $ 通过栅格的平流项模拟实现。3.2 栅格的平流项模拟
以中心栅格为单元，计算时间片上邻域栅格流体流向中心栅格的水量，加上中心栅格剩余的水量作为时间片模拟完成后中心栅格的水量，以模拟式（5）中的
$ In_r^n $ 和$ Out_r^n $ 。设中心栅格c的8邻域栅格为b，其水流速度矢量为
$ \left( {v_b^x,v_b^y} \right) $ ，则邻域栅格流体流向中心栅格的水量比例$S _{b \to c}$ 为：$$ \begin{split} &S _{b \to c} = \\ &\max \left( {1 - \left| {\left. {\varDelta x + v_b^x/{v_{\max }}} \right|,0} \right.} \right) \times \max \left( {1 - \left| {\left. {\varDelta y + v_b^y/{v_{\max }}} \right|,0} \right.} \right) \end{split} $$ (7) 式中：
$\left( {\varDelta x,\varDelta y} \right)$ ——单位像素栅格之间所产生的坐标偏差；$ {v_{\max }} $ ——水流流动时在栅格相等时间片上最大速度 限制。中心栅格剩余水量比例
$S _c$ 为：$$ S _c = \max \left( {1 - \left| {\left. {v_c^x/{v_{\max }}} \right|,0} \right.} \right) \times \max \left( {1 - \left| {\left. {v_c^y/{v_{\max }}} \right|,0} \right.} \right) $$ (8) 式7和式8中，
$ v_{b}^{t} $ 和$ v_{t}^{y} $ 的速度标量不能大于${v_{\max }} $ ，否则将导致水量失衡。当模拟的水速标量大于${v_{\max }} $ 时，需强制设置为$ {v_{\max }} $ 。因此，${v_{\max }} $ 为模型模拟的最大速度限制，决定了时间片模拟长度，即：$$ \varDelta t = C/{v_{\max }} $$ (9) 式中：C——栅格分辨率；
$ \varDelta t $ ——时间片长度。基于
$S _{b \to c}$ 和$S _c$ ，时间片n+1的中心栅格水深为：$$ W_c^{n + 1} = \left( {\sum\limits_{b = 1}^8 {W_b^nS _{b \to c}^n} } \right) + W_c^nS _c^n $$ (10) 式中：
$ W_b^n $ ——中心栅格邻域时间片的初始水深；$ S _{b \to c}^n $ ——时间片上邻域栅格b在流体运动时进入到 中心栅格c中的水量比例；$ S_c^n $ ——时间片上中心栅格剩余的水量比例。设
$ V_c^n $ 和$ V_c^{n + 1} $ 分别表示中心栅格时间片上的初始水速和终止水速。依据动量守恒定律，终止水速为：$$ V_c^{n + 1} = \frac{{\displaystyle\sum\limits_{b = 1}^8 {\dfrac{{{{\left( {V_{b \to c}^n} \right)}^2}\varDelta t}}{{\varDelta L}}W_b^n + V_c^n\left( {W_c^n - Out_c^n} \right)} }}{{W_c^{n + 1}}} $$ (11) 3.3 栅格的压力项和外力项模拟
栅格的水动力模拟包括压力项和外力项模拟。地表水在流动过程中受重力、压力等的影响，产生相应的加速度，用
$ \varDelta {H_b} $ 来代表单元之间水体流动的高程差，邻域栅格对中心栅格所产生的速度矢量$ \left({v}_{b}^{x}，\varDelta {v}_{b}^{y}\right) $ 为：$$ \left\{ \begin{gathered} \varDelta v_b^x = \sum _{\left( {j,i} \right)}{\left( { - j \cdot \alpha \cdot \varDelta {H_b}} \right)} \\ \varDelta v_b^y = \sum _{\left( {j,i} \right)}{\left( { - i \cdot \alpha \cdot \varDelta {H_b}} \right)} \\ \end{gathered} \right. $$ (12) 式中：j、i——水流流动模拟过程中邻域栅格与中心栅格的列序号和行序号差；
$ \alpha $ ——与加速度、水的密度、摩擦力等相关的率定正数常量。水体在地表流动的过程中时会受地表摩擦力影响而使水流速度发生纵向的变化，当实际水深低于给定的阈值时，设定一个衰减系数以减小模拟水速，衰减系数
$ \varepsilon $ 近似定义为：$$ \varepsilon = \left\{ \begin{split} & \sigma \cdot {\left( {\frac{{{W_c} - {d_{\min }}}}{{{d_{\max }} - {d_{\min }}}}} \right)},&{W_c} \leqslant {d_{\max }} \\ & {\sigma }\quad\quad\quad\quad\quad\;\;\;,&{W_c} > {d_{\max }} \\ \end{split} \right. $$ (13) 式中：
$ {d_{\min }} $ 、$ {d_{\max }} $ —水体流动所受摩擦力的下界深度和 上界深度/m；$ \sigma $ ——流体的摩擦比例系数。3.4 冰川融雪模拟
林芝市区域内有大面积冰川分布，其高山地区为典型的海洋性冰川，同时受季风影响，所以积雪消融速度也相对较快，冰川消融型泥石流是林芝市泥石流灾害的主要类型之一。冰川消融所带来的水分也为周围流域沟谷形成泥石流灾害提供了一定的基础。本文将冰川消融加入径流模型模拟中，作为径流模拟初始水源部分。前人大都采用能量平衡[23-24]的方法对冰雪消融进行模拟。能量平衡的计算复杂并且需要大量资料，由于西藏林芝地势地形特殊、资料稀缺，能量平衡算法在当地难以推广。度日模型能够在流域尺度上基于气温模拟出与能量平衡模型相似的结果，并且相比能量平衡模型数据更易获取、计算更为简单。度日模型基本计算公式如下：
$$ M = \left\{ \begin{split} &MF(T - {T_0}),&{T > {T_0}} \\ &0\quad\quad\quad\quad\;\;,&{T \leqslant {T_0}} \\ \end{split} \right. $$ (14) 式中：
$ M $ ——日消融量/mm；$ T $ ——日平均气温/°C；$ {T_0} $ ——冰川消融温度阈值/°C；$ MF $ ——度日因子/（mm·°C−1·d−1）。将气温数据插值后代入式（14）计算得到冰川所在栅格像元的日消融量，将融冰量加入了净雨量中[25-26]，再进行产汇流计算。
3.5 径流汇流模型的精度评价与水深分级
选择2015年8月19日00时为起始时间，模拟时长18天，将气象站的降水和气温数据代入径流汇流模型，取
$ \nu_{\max } $ =2，$ \alpha $ =0.24，$ \sigma $ =0.95，以栅格空间分辨率250 m模拟林芝市的径流汇流过程。图4为2015年9月6日23时的林芝市径流汇流模拟结果。从图中可见，本文模型可获取流域各位置上的水深，并且水深较大区域为沟谷和河流。![]()
图 4 2015年9月6日23时径流汇流模拟结果Figure 4. Results of runoff simulation at 23:00 on September 6, 2015为验证模拟的合理性，在逐小时的径流汇流模拟栅格数据上，采集雅鲁藏布江河段下游水文站（图4）的模拟水深，将实际观测水深与模拟水深进行对比。图5为模拟水深与实际水深对比曲线。从图中可见，模拟水深曲线与实际水深曲线趋势一致。将实际水深数据与模拟水深数据进行pearson相关性分析，其相关性为0.838，为强正相关，表明模型能模拟出水深的变化趋势。
![]()
图 5 测量水位点模拟水深与实际水深对比Figure 5. Comparison between simulated water depth and actual water depth at measured water level sites虽然模型模拟的水深变化趋势接近，但模拟水深的变化值与实际水深变化值不在相同量级上。这是因为，本文的径流汇流在干涸河道条件上进行的，未考虑河道初始水位；同时，模型输出的水深是以栅格像元为单位的，相当于将像元内的水平铺于250×250 m2的矩形栅格单元上。综上，以径流汇流模型为基础的泥石流预警模型，需要利用实测泥石流灾害点和模拟水深进行相关性分析，确定模拟水深的权重。
将林芝市泥石流灾害点与径流汇流模拟水深进行叠加分析，发现灾害点模拟水深均大于0.01 m，水深在0.01～0.05 m的灾害点占比8%，水深在0.05～0.1 m的灾害点占比12%，水深在0.1～0.3 m的灾害点占比42%，水深在0.3 m以上的灾害点占比38%。将水深按表3分级，作为泥石流预警模型的输入。
表 3 径流水深分级Table 3. Classification of runoff depth级别 水深区间/m 1 (0, 0.01] 2 (0.01, 0.05] 3 (0.05, 0.1] 4 (0.1, 0.3] 5 >0.3 4. 泥石流预警模型
4.1 泥石流灾害预警方法
通常而言，泥石流从物源区启动。泥石流是在物源区以下均存在。从防灾减灾工程而言，更关注人类活动的区域。该区域内，容易引发经济财产损失和威胁生命安全。林芝市主要城镇、村庄、农田、草场、公路大多分布在江河两岸的阶地与滩地上或谷地中。因此，林芝市历史泥石流灾害点一般标注在沟谷、沟口及附近。为进一步提高灾害预测的空间针对性，将流域出口点作为沟口，以500 m为缓冲距离，对沟口距离进行反距离归一化计算：
$$ L = \left\{ \begin{split} &(500 - d)/500,&{d \leqslant 500} \\ & 0\quad\quad\quad\quad\quad\;,&{d > 500} \end{split} \right. $$ (15) 利用式15，产生林芝市沟口距离归一化栅格图，并将其与历史128个泥石流灾害点进行叠加统计。结果表明，沟口距离归一化值0.8～0.6范围内灾害点共18个，占比14.29%，在0.8～1范围内灾害点共79个，占比63.7%（图6）。
![]()
图 6 林芝市泥石流历史灾害点河道距离Figure 6. The distance of river course of historical disaster site of debris flow in Nyingchi City综合径流汇流模拟水深分级
$ D $ ，泥石流灾害危险性概率$ P $ ，沟口距离归一化因子$ L $ ，利用二元逻辑回归的方法计算各个因素对泥石流灾害的计算权重，如下公式：$$ Y = \left( {aD + bP} \right) \cdot L $$ (16) 式中：
$ Y $ ——林芝市泥石流灾害预警值；$ D $ ——水深分级；$ P $ ——泥石流灾害危险性概率；$ a、b $ ——权重。4.2 基于逻辑回归的泥石流灾害预警模型
选取2000—2015年的历年泥石流灾害点数据，参与逻辑回归运算，以泥石流灾害是否发生作为因变量，径流汇流模拟水深、泥石流灾害背景值作为自变量，二元逻辑回归的权重系数见表4：
表 4 逻辑回归方程表Table 4. Logistic regression equation tableB S.E. Wald df Sig. exp（B） P 0.790 0.257 9.486 1 0.002 2.204 D 0.834 0.297 7.910 1 0.005 2.303 常量 −5.330 1.432 13.857 1 0.000 0.005 注：B为逻辑回归系数；S.E.为标准误差；Wald为卡方值；df为自由度；Sig.为显著性。 对泥石流灾害指标进行逻辑回归后，其所有因子显著性都小于0.05，基本可判定该逻辑回归模型具有显著性，逻辑回归模型通过检验。将B值归一化，作为预警因子的权重：
$$ Y = \left( {{\text{0}}{\text{.52}}D + {\text{0}}{\text{.48P}}} \right) \times L $$ (17) 以2011—2016年33次泥石流灾害为研究对象，利用本文模型计算得到Y值分布图。首先，设预警阈值为a，获取≥a的区域作为预警区域，将预警区域面积除以区域总面积，得到单次预警区面积占比，统计33次灾害不同阈值下的平均预警区面积平均占比；然后，与实际灾害点进行叠加，统计落入预警区的灾害点个数，获取不同阈值下预警区域内的灾害点占比。统计结果如图7所示。
![]()
图 7 不同预警阈值下的预警区灾害点比例和平均预警区面积平均占比Figure 7. Proportion of disaster sites in early warning area and average proportion of average early warning area under different early warning thresholds本文模型融入径流汇流算法后，预警区绝大部分分布在河流和沟谷地区。从图7可见，预警阈值为0.15时，平均预警区面积平均占比为0.05%，随着阈值a的增加，预警区面积占比不断减少，这反应了模型较强的空间针对性。为增加落入预警区的灾害点数量，期望较大的预警阈值，但较大的预警阈值将导致预警区面积太大使预警结果失去实际意义。研究发现，当预警阈值为0.15时，河流和沟谷区域基本上为预警区域，模型的灾害风险辨别能力较低。当阈值a在0.55附近时，图7中的预警区面积占比突然下降至0.005%，但预警精度仍保持在60%左右。因此，综合预警区面积占比和预报精度，将0.55作为模型的预警阈值。
4.3 泥石流灾害预警精度评价
以0.55为预警阈值，以2011—2020年历史灾害为结果验证数据，得到验证结果如表5。
表 5 泥石流灾害预警验证结果Table 5. Validation results of debris flow disaster early warning序号 灾害日期 灾害点数/个 预警区内灾害点数/个 1 2011年7月14日 2 1 2 2013年7月6日 1 1 3 2014年4月16日 1 0 4 2014年7月18日 1 0 5 2015年8月6日 5 1 6 2015年8月19日 12 7 7 2015年8月20日 2 2 8 2015年8月21日 3 3 9 2016年4月25日 1 1 10 2016年7月26日 3 3 11 2016年7月27日 2 1 12 2017年8月3日 1 0 13 2017年7月21日 1 1 14 2018年9月9日 1 1 15 2019年7月7日 2 1 16 2020年7月10日 2 2 17 2020年7月11日 3 2 18 2020年8月15日 2 2 合计 45 29 从表5可见，案例对应的历史灾害点数据共45个，预警区内灾害点数为29个，占比64.4%。2016年后，历史灾害点数据共18个，预警区内灾害点数为14个，占比77.8%。2016年以前未预测到灾害点大部分原因是气象站较少，降雨数据缺失，导致部分地区径流汇流模拟结果欠佳即预警效果降低。自2016年以后由于开始逐渐新增气象站点，预警模型精度也在一定程度上得到了提高。
图8为林芝市2020年7月11日和8月15日的泥石流灾害预警图。7月11日、8月15日发生泥石流灾害均为1处，灾害点发生处预报为红色，表明本文的泥石流预警在具有相对高的空间针对性的同时，能很好地反应泥石流发生的风险。
![]()
图 8 林芝市2020年7月11日、8月15日泥石流灾害预警图Figure 8. Debris flow disaster warning map of Nyinginzhi on July 11、August 15, 20205. 结论与讨论
5.1 结论
为解决传统基于降水量的泥石流预警模型无法预测泥石流灾害发生部位，提出了一种基于栅格径流汇流模拟的预警模型。该模型将泥石流致灾因子分成了2类，即背景因子和激发因子。选取裸岩率、流域面积、纵比降、河流、道路、断层密度、土地利用、土壤类型、隐患点密度、沟谷密度、年降水量等，利用逻辑回归计算得到林芝市泥石流灾害概率，作为模型的背景因子；以径流水深作为模型的激发因子，构建了栅格模式下的径流汇流模型。通过栅格时间片上的平流项、压力项和外力项模拟，计算得到流域各位置的模拟水深。研究发现，模型模拟水深与实测水深虽然不在同一个量级，但较好地反应了水深变化趋势。因此，将模拟水深进行分级，并综合沟口距离和泥石流灾害概率，构建林芝市泥石流灾害预警模型。利用2011—2020年18次历史灾害点对模型进行验证，预报等级在黄色以上灾害点为29次，占比64.4%。
5.2 讨论
在本文的径流汇流模拟中，由于缺少流域前期水情资料，无法构建流域前期基础流场。同时，径流汇流模拟在250 m的空间分辨率上进行的，栅格单元的水被分摊在250×250 m2的栅格单元上，导致模拟水深与实际观测水深不在同一量级上。通过实际灾害数据进行分析，发现泥石流灾害均发生在模拟水深0.01 m以上的区域，其中有80%的灾害发生在模型模拟水深0.1 m以上区域，说明径流模拟模型在一定程度上能反应了林芝地区泥石流灾害发生的概率。为合理融入径流汇流的水深，将模拟水深依据历史灾害进行分级，并将分级值导入模型，采用逻辑回归法相对客观地确定模拟水深的权重因子。实验证明，该方法在保证一定的泥石流预警精度的同时，避免了水深模拟误差对预警精度的直接影响。
另外，数据精度决定模型精度。气象站数量、面雨量插值方法等均对模型精度有直接的影响。2015年之前，由于气象站数量较少，面雨量精度较低，模型预警精度较低；2016年以后，由于林芝市站点的增加，预警精度达到77.8%。径流汇流模型需要输入区域面雨量，由于无法直接获取高精度的网格降水，采用空间插值方法近似表示区域的面雨量，加之地形数据本身误差，均导致了模型精度难以进一步提高。
从模型构建方法上来看，虽然本文计算出来的预警值能反应区域泥石流风险状况，但阈值法采用了较为粗略的分析方法，可能给预警结果产生不确定性的影响。另外，由于数据资料有限，2020年以来仅有3次泥石流灾害资料，无法说明模型的精度，只能将验证时段加长，采用2011—2020年的历史灾害进行精度验证。需要进一步增加验证资料，并不断完善模型，以进一步提高预警精度。
-
![]()
图 1 林芝市泥石流灾害点、隐患点分布图
Figure 1. Distribution of debris flow disaster sites and hidden danger points in Nyingchi City
![]()
图 4 2015年9月6日23时径流汇流模拟结果
Figure 4. Results of runoff simulation at 23:00 on September 6, 2015
![]()
图 5 测量水位点模拟水深与实际水深对比
Figure 5. Comparison between simulated water depth and actual water depth at measured water level sites
![]()
图 6 林芝市泥石流历史灾害点河道距离
Figure 6. The distance of river course of historical disaster site of debris flow in Nyingchi City
![]()
图 7 不同预警阈值下的预警区灾害点比例和平均预警区面积平均占比
Figure 7. Proportion of disaster sites in early warning area and average proportion of average early warning area under different early warning thresholds
![]()
图 8 林芝市2020年7月11日、8月15日泥石流灾害预警图
Figure 8. Debris flow disaster warning map of Nyinginzhi on July 11、August 15, 2020
表 1 逻辑回归分析结果
Table 1 Results of logistic regression analysis
指标因子 B S.E Wals df Sig. exp（B） 裸岩信息 0.925 0.456 4.119 1 0.042 2.523 流域面积 0.880 0.394 5.000 1 0.025 2.412 沟床纵比降 1.042 0.186 31.348 1 0 2.834 河流 1.000 0.122 66.763 1 0 2.717 道路 0.879 0.105 69.489 1 0 2.409 断层密度 0.995 0.442 5.059 1 0.024 2.704 土地利用 0.541 0.245 4.888 1 0.027 1.718 土壤类型 0.603 0.112 18.239 1 0 1.828 隐患点密度 0.896 0.297 63.593 1 0 2.449 沟谷密度 1.295 0.211 19.043 1 0 3.650 年降水量 1.064 0.281 25.459 1 0 2.897 常量 −1.485 0.281 27.852 1 0 0.226 注：B为逻辑回归系数；S.E.为标准误差；Wals为卡方值统计量；df为自由度；Sig.为显著性。 
下载: 导出CSV
表 2 泥石流灾害的概率占比
Table 2 Probability proportion of debris flow disaster
级别 取值区间 灾害点数量 百分比/% 一级 [0, 0.2) 8 6.25 二级 [0.2, 0.4) 5 3.90 三级 [0.4, 0.6) 2 1.56 四级 [0.6, 0.8) 9 7.03 五级 [0.8, 1) 104 81.25
下载: 导出CSV
表 3 径流水深分级
Table 3 Classification of runoff depth
级别 水深区间/m 1 (0, 0.01] 2 (0.01, 0.05] 3 (0.05, 0.1] 4 (0.1, 0.3] 5 >0.3
下载: 导出CSV
表 4 逻辑回归方程表
Table 4 Logistic regression equation table
B S.E. Wald df Sig. exp（B） P 0.790 0.257 9.486 1 0.002 2.204 D 0.834 0.297 7.910 1 0.005 2.303 常量 −5.330 1.432 13.857 1 0.000 0.005 注：B为逻辑回归系数；S.E.为标准误差；Wald为卡方值；df为自由度；Sig.为显著性。
下载: 导出CSV
表 5 泥石流灾害预警验证结果
Table 5 Validation results of debris flow disaster early warning
序号 灾害日期 灾害点数/个 预警区内灾害点数/个 1 2011年7月14日 2 1 2 2013年7月6日 1 1 3 2014年4月16日 1 0 4 2014年7月18日 1 0 5 2015年8月6日 5 1 6 2015年8月19日 12 7 7 2015年8月20日 2 2 8 2015年8月21日 3 3 9 2016年4月25日 1 1 10 2016年7月26日 3 3 11 2016年7月27日 2 1 12 2017年8月3日 1 0 13 2017年7月21日 1 1 14 2018年9月9日 1 1 15 2019年7月7日 2 1 16 2020年7月10日 2 2 17 2020年7月11日 3 2 18 2020年8月15日 2 2 合计 45 29
下载: 导出CSV
-
[1] 屈永平, 唐川, 刘洋, 等. 西藏林芝地区冰川降雨型泥石流调查分析[J]. 岩石力学与工程学报, 2015, 34（增刊2）: 4013 − 4022 QU Yongping, TANG Chuan, LIU Yang, et al. Investigation and analysis of glacier debris flow in Nyingchi area, Tibet[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34（Sup 2）: 4013 − 4022. （in Chinese with English abstract）
[2] 陈宫燕,普布桑姆,次仁旺姆,等. 基于Logistic回归方法的林芝市山洪地质灾害预警研究[J]. 中国农学通报,2019,35(23):124 − 130. [CHEN Gongyan,PU Busangmu,CI Renwangmu,et al. Early warning of mountain flood geological disasters in Nyingchi based on logistic regression method[J]. Chinese Agricultural Science Bulletin,2019,35(23):124 − 130. (in Chinese with English abstract) DOI: 10.11924/j.issn.1000-6850.casb18040021 [3] 陈宫燕,普布桑姆,次仁,等. 西藏林芝降水引发的山洪地质灾害分布特征[J]. 中国地质灾害与防治学报,2018,29(2):100 − 103. [CHEN Gongyan,PU Busangmu,CIREN,et al. Distribution characteristics of mountain flood and geological disaster caused by precipitation of Nyingchi in Tibet[J]. The Chinese Journal of Geological Hazard and Control,2018,29(2):100 − 103. (in Chinese with English abstract) [4] SALINAS-JASSO J A,VELASCO-TAPIA F,DE LEÓN I N,et al. Estimation of rainfall thresholds for shallow landslides in the Sierra Madre Oriental,northeastern Mexico[J]. Journal of Mountain Science,2020,17(7):1565 − 1580. DOI: 10.1007/s11629-020-6050-2
[5] YAN Y,ZHANG Y,HU W,et al. A multiobjective evolutionary optimization method based critical rainfall thresholds for debris flows initiation[J]. Journal of Mountain Science,2020,17(8):1860 − 1873. DOI: 10.1007/s11629-019-5812-1
[6] 马煜,李彩侠. 地形因子对四川龙溪河流域泥石流发生的影响[J]. 中国地质灾害与防治学报,2018,29(2):16 − 21. [MA Yu,LI Caixia. Effect of topographic factor on the formation of debris flows in Longxi River Area,Sichuan Province[J]. The Chinese Journal of Geological Hazard and Control,2018,29(2):16 − 21. (in Chinese with English abstract) [7] 熊江,唐川,陈明. 泥石流早期识别与监测预警研究进展探讨[J]. 自然灾害学报,2021,30(1):165 − 173. [XIONG Jiang,TANG Chuan,CHEN Ming. Discussion on the research progress of early identification,monitoring and early warning of debris flow[J]. Journal of Natural Disasters,2021,30(1):165 − 173. (in Chinese with English abstract) [8] 王英杰. 岷江上游泥石流预警报模型构建及其应用[D]. 绵阳: 西南科技大学, 2020 WANG Yingjie. Construction and application of debris flow early-warning and forecast model in the upper Min River[D]. Mianyang: Southwest University of Science and Technology, 2020. （in Chinese with English abstract）
[9] PONZIANI M,POGLIOTTI P,STEVENIN H,et al. Debris-flow Indicator for an early warning system in the Aosta valley region[J]. Natural Hazards,2020,104(2):1819 − 1839. DOI: 10.1007/s11069-020-04249-5
[10] 黄健,胡卸文,金涛,等. 四川西昌“3·30”火烧区响水沟火后泥石流成灾机理[J]. 中国地质灾害与防治学报,2022,33(3):15 − 22. [HUANG Jian,HU Xiewen,JIN Tao,et al. Mechanism of the post-fire debris flow of the Xiangshui gully in “3·30” fire area of Xichang, Sichuan Province[J]. The Chinese Journal of Geological Hazard and Control,2022,33(3):15 − 22. (in Chinese with English abstract) [11] 包红军,王凯,张少杰,等. 耦合分布式水文模型的泥石流物理模型预报试验[J]. 暴雨灾害,2018,37(4):303 − 310. [BAO Hongjun,WANG Kai,ZHANG Shaojie,et al. Debris flow forecasting test based on a physical model coupling with a distributed hydrological model[J]. Torrential Rain and Disasters,2018,37(4):303 − 310. (in Chinese with English abstract) [12] 王治华,郭兆成,杜明亮,等. 基于数字滑坡技术的暴雨滑坡、泥石流预警、监测模型研究[J]. 地学前缘,2011,18(5):303 − 309. [WANG Zhihua,GUO Zhaocheng,DU Mingliang,et al. Model study of monitoring and early warning of rainstorm induced landslide and debris flow based on digital landslide technology[J]. Earth Science Frontiers,2011,18(5):303 − 309. (in Chinese with English abstract) [13] 王海芝. 北京地区暴雨泥石流预警阈值研究[J]. 第四纪研究,2020,40(5):1371 − 1380. [WANG Haizhi. Early warning thresholds of the rainfall-induced debris flows in Beijing[J]. Quaternary Sciences,2020,40(5):1371 − 1380. (in Chinese with English abstract) DOI: 10.11928/j.issn.1001-7410.2020.05.24 [14] BUGAETS A N,GARTSMAN B I,TERESHKINA A A,et al. Using the SWAT model for studying the hydrological regime of a small river basin (the komarovka river,primorsky Krai)[J]. Russian Meteorology and Hydrology,2018,43(5):323 − 331. DOI: 10.3103/S1068373918050060
[15] THOMPSON J R,SØRENSON H R,GAVIN H,et al. Application of the coupled MIKE SHE/MIKE 11 modelling system to a lowland wet grassland in southeast England[J]. Journal of Hydrology,2004,293(1/2/3/4):151 − 179.
[16] 王济,蔡雄飞,雷丽,等. 不同裸岩率下我国西南喀斯特山区土壤侵蚀的室内模拟[J]. 中国岩溶,2010,29(1):1 − 5. [WANG Ji,CAI Xiongfei,LEI Li,et al. Laboratory simulation on soil erosion under different bedrock outcrop rate in Southwest Karst area, China[J]. Carsologica Sinica,2010,29(1):1 − 5. (in Chinese with English abstract) [17] LIU W, HE S. Comprehensive modelling of runoff-generated debris flow from formation to propagation in a catchment. Landslides, 2020, 17（4）: 1529 − 1544.
[18] 张书豪,吴光,张乔,等. 基于子流域特征的泥石流易发性评价[J]. 水文地质工程地质,2018,45(2):142 − 149. [ZHANG Shuhao,WU Guang,ZHANG Qiao,et al. Debris-flow susceptibility assessment using the characteristic factors of a catchment[J]. Hydrogeology & Engineering Geology,2018,45(2):142 − 149. (in Chinese with English abstract) [19] 张以晨,秦胜伍,翟健健,等. 基于信息量的长白山地区泥石流易发性评价[J]. 水文地质工程地质,2018,45(2):150 − 158. [ZHANG Yichen,QIN Shengwu,ZHAI Jianjian,et al. Susceptibility assessment of debris flow based on GIS and weight information for the Changbai Mountain area[J]. Hydrogeology & Engineering Geology,2018,45(2):150 − 158. (in Chinese with English abstract) [20] 张盼盼,胡远满,肖笃宁,等. 一种基于多光谱遥感影像的喀斯特地区裸岩率的计算方法初探[J]. 遥感技术与应用,2010,25(4):510 − 514. [ZHANG Panpan,HU Yuanman,XIAO Duning,et al. A method of the percentage of bare rock calculation in Karst areas based on multi-spectrum remote sensing image[J]. Remote Sensing Technology and Application,2010,25(4):510 − 514. (in Chinese with English abstract) DOI: 10.11873/j.issn.1004-0323.2010.4.510 [21] 陈冠,孟兴民,郭鹏,等. 白龙江流域基于GIS与信息量模型的滑坡危险性等级区划[J]. 兰州大学学报(自然科学版),2011,47(6):1 − 6. [CHEN Guan,MENG Xingmin,GUO Peng,et al. Landslide susceptibility mapping based on GIS and information value model in Bailong River Basin[J]. Journal of Lanzhou University (Natural Sciences),2011,47(6):1 − 6. (in Chinese with English abstract) [22] 王卫东,钟晟. 基于GIS的Logistic回归模型在地质灾害危险性区划中的应用[J]. 工程勘察,2009,37(11):5 − 10. [WANG Weidong,ZHONG Sheng. The application of GIS-based logistic regression for geological hazards zonation[J]. Geotechnical Investigation & Surveying,2009,37(11):5 − 10. (in Chinese with English abstract) [23] HOCK R. Glacier melt:A review of processes and their modelling[J]. Progress in Physical Geography:Earth and Environment,2005,29(3):362 − 391. DOI: 10.1191/0309133305pp453ra
[24] 张勇,刘时银. 度日模型在冰川与积雪研究中的应用进展[J]. 冰川冻土,2006,28(1):101 − 107. [ZHANG Yong,LIU Shiyin. Progress of the application of degree-day model to study glaciers and snow cover[J]. Journal of Glaciology and Geocryology,2006,28(1):101 − 107. (in Chinese with English abstract) DOI: 10.3969/j.issn.1000-0240.2006.01.014 [25] 闫玉娜,车涛,李弘毅,等. 使用积雪遥感面积数据改善山区春季融雪径流模拟精度[J]. 冰川冻土,2016,38(1):211 − 221. [YAN Yuna,CHE Tao,LI Hongyi,et al. Using snow remote sensing data to improve the simulation accuracy of spring snowmelt runoff:Take Babao River Basin as an example[J]. Journal of Glaciology and Geocryology,2016,38(1):211 − 221. (in Chinese with English abstract) [26] OMANI N,SRINIVASAN R,SMITH P K,et al. Glacier mass balance simulation using SWAT distributed snow algorithm[J]. Hydrological Sciences Journal,2017,62(4):546 − 560. DOI: 10.1080/02626667.2016.1162907
-
期刊类型引用(1)
1. 孟迎. 考虑降雨量数值变化的山洪灾害动态预警方法研究. 水资源开发与管理. 2024(03): 73-76+56 . 
百度学术
其他类型引用(2)
计量
- 文章访问数: 1977
- HTML全文浏览量: 1667
- PDF下载量: 240
- 被引次数: 3
邮件订阅
RSS