基于水动力模型的泥石流风险分析

张弘强， 李文娟， 李大鸣， 唐智超， 汤稀博， 唐颖

(1.国网湖南省电力有限公司电力科学研究院， 长沙 410007；

2.高效清洁发电技术湖南省重点实验室， 长沙 410007；

3.湖南省湘电试验研究院有限公司， 长沙 410007；

4.天津大学水利工程仿真与安全国家重点实验室， 天津 300072；

5.国网江西省电力有限公司萍乡市湘东区供电分公司， 萍乡 337016；

6.河海大学能源与电气学院， 南京 211100；

7.长沙理工大学水利工程学院， 长沙 410114)

泥石流是一种破坏力强的常见灾害，因此逐渐引起了中外学者的关注。对区域泥石流危险度进行评估，对泥石流风险进行分析，对于减轻泥石流灾害的损失具有重要意义[1-2]。Mizuyama等[3]提出了用于计算由泥石流形成的沟床的比降变化的一维运动方程，以及用于计算泥石流影响范围的二维运动方程。O’Brien等[4]将适用于黏性泥石流的Bing模型以及适用于稀性泥石流的Bagnold模型做了合并处理，提出了一种通用的泥石流运动模型，即稀释塑性模型。De Joode等[5]构建了一种增加了暴雨影响的泥石流灾害模型，并将其与地理信息软件相结合。Martinez等[6]提出了一种近似三维的数值模拟方法，用于模拟碎石较多的泥石流。Zou等[7]提出一种基于水文响应单元定量评价泥石流区域风险的方法。国内学者从20世纪末开始对泥石流危险区域进行划分，并得到了一系列相当具有价值的成果。谭炳炎[8]对泥石流沟进行了实地调查和统计，最先对泥石流沟的严重程度进行了数量化综合评价。刘希林等[9]通过分析和调查，确定泥石流危险度的主要危险因子和次要危险因子，根据关联次序确定它们的权重，以此提出了泥石流危险度的定量计算公式，极大地推动了泥石流量化评价的发展。李大鸣等[10]利用模糊赋值法计算各个泥石流危险因子的权重，建立了天津蓟县泥石流区域预报模型，确定了天津蓟县泥石流区域危险度。陈飞飞等[11]对泥石流危险评价问题进行探讨，指出了目前评价方法中存在的问题，并探讨了相应的解决方法。王高峰等[12]利用MATLAB多元非线性统计方法建立白龙江流域泥石流一次最大冲出量预测模型。李丽敏等[13]提出基于布谷鸟算法优化后的轻量梯度提升机的泥石流预测。当前的泥石流评价主要集中在泥石流的静态评价，在泥石流动态评价方面涉及较少[11]。

为构建泥石流动态评价机制，现以水动力模型为研究依据，以降水、来流等自然条件为泥石流发生的推动因子，加之地形、地质、植被、人口等影响因子，构建泥石流风险分析模型。为准确判断泥石流灾害危险区域分布，为减轻泥石流损害提供指导性技术，现选取广西5个市的行政区域作为研究区域，通过构建水动力数学模型，以降水、地形、地质、植被、人口等作为泥石流的影响因子，采用河网分级编码方法概化河网，采用粗糙集理论计算权重，从而进行泥石流风险分析，并将计算模拟结果与广西地质灾害分布统计和2011年全州县泥石流灾害事件进行对比，验证模型的合理性和可靠性。

1.1 控制方程

若假设在水深方向的动水压强分布规律与静水压强相一致，则平面二维连续方程[14]为

(1)

式(1)中：t为时间；
h为水深；
u、v分别为x、y方向上的速度平均值；
q为源汇项。

采用莱布尼兹公式当中的分解形式，并对其积分项作合并，从而得到二维浅水运动方程为

(2)

式(2)中：n为糙率；
g为重力加速度。

1.2 控制方程离散求解

采用有限体积法对式(2)作离散，对于时间层面上的计算方式，采取水位与流量交错的计算方式。对于控制方程的通用形式，可表示为

(3)

式(3)中：ρ为水的密度；
φ为广义变量；
Γ为φ所对应的广义扩散系数；
S为φ所对应的广义源项。

至于时间段Δt与控制体积ΔV，对控制方程的通用式(3)作积分可得

(4)

式中：div为求矢量的散度；
grad为求矢量的梯度。

对于式(4)在空间与时间上的积分处理，集中表现在以下4个方面。

(1)瞬态项积分。

(5)

(2)源项积分。

(6)

式(6)中：SC对应节点C上广义变量所对应的广义源项的大小；
SP对应节点P上广义变量所对应的广义源项的大小；
φP为对应节点P上广义变量的大小。

(3)对流项积分。

(7)

式(7)中：u为速度；
(ρuφA)e、(ρuφA)w、(ρvφA)n和(ρvφA)s分别为ρuφA在所对应节点e、w、n和s的大小；
(ρu)e、(ρu)w分别为ρu在所对应节点e、w的大小；
(ρv)w、(ρv)s为ρv在所对应节点w、s的大小；
Ae、Aw、An和As分别为A在所对应节点e、w、n和s的大小；
φe、φw、φn和φs分别为φ在所对应节点e、w、n和s的大小。

(4)扩散项积分。

(8)

在对连续方程、动量方程等作离散求解时，只需要把相对应的物理量作为φ代入式(4)～式(8)，即可得到相对应的离散方程。对于连续微分方程，其对任何多边形网格的显式离散表达式为

(9)

式(9)中：Hi为i号网格单元的模化水深；
T为所处时刻；
Ai为i号网格单元所占面积大小；
Lik为i号网格、k号通道的长度；
Qik为i号网格、k号通道的单宽流量；
qi为i号网格的抽排水量。

1.3 边界条件确定

对于模型边界条件的确定，按照以下3个方面进行选取。

(1)研究区域边界条件。以桂林、贺州、梧州、来宾、柳州5个市作为研究区域，以外轮廓线作为模型区域边界条件。

(2)降水边界条件。将开始时刻陆域网格内的水深大小设为零，采用曲面差值将测量站所采集到的点雨量数据转化为面雨量数据，以此构成模拟计算泥石流的降水边界条件。

(3)下垫面边界条件。将广西的地形、地质、植被、人口分布等纳入计算范围，生成泥石流模型的下垫面边界条件。

1.4 临界条件确定

从泥石流的产生条件与影响因子出发，基于降水强度与入渗率的比较，分析泥石流的起动过程，即：①降水强度小于入渗率，雨水入渗，无积水；
②浅层土体开始饱和，降水强度等于或略大于入渗率，积水导致浅层径流，开始冲蚀松散土体；
③径流冲刷松散土体，地下渗流导致土体间的吸引力减小、抗剪强度减小，从而导致剪切破坏，加之水流的不断作用，最终造成泥石流灾害的发生。

根据团队之前研究成果，给出泥石流发生的临界条件[15-16]为

hc=[(ρsatgHcosθ+ρgH)(tanφ-tanθ)-ρsgHtanφ+τ1+τ2]/[ρmgcosθ(2tanθ-tanφ)+ρgtanφ]

(10)

式(10)中：hc为临界径流深度；
H为破坏面入渗高度；
ρsat为松散土体的饱和密度大小；
ρm为松散固体物质的密度大小；
ρs为坡面松散固体物质密度；
ρ为水的密度；
τ1为土体性质的抗冲蚀能力的概化；
τ2为植被覆盖的抗冲蚀能力的概化；
θ为斜坡坡度；
φ为内摩擦角。

1.5 泥石流数学模型构建

图1 广西数字高程模型Fig.1 Digital elevation model of Guangxi

广西境内山地较多，地势较为陡峭，高低落差较大，加之广西地区雨量充沛，降水分布一般呈现出山区丘陵地带多、平原河谷地带少的规律。尤其在夏季，迎风坡面降水多、背风坡面降水少的特点，无疑为泥石流的发生创造了先决条件。从广西数字高程模型图(图1)中提取整个研究区域的边界，并将研究区域内的35个市县的地理信息数字化，得到研究区域的行政分区(图2)。

图2 模型行政区域划分Fig.2 Administrative division of model

1.5.1 网格划分

对于网格划分，采用不规则、无结构网格划分技术，对研究区域进行自动剖分、逐步加密与分块优化。为了能够确定研究区域占矩形范围的大小，以边长为0.01°的曲面正方形网格对矩形范围内的网格进行自动划分，同时分别对单元、节点与通道进行编号。研究区域网格单元大小为2 km×2 km，模型共生成网格22 900个、节点23 768个、通道46 667个，局部示意如图3所示。

1.5.2 地形处理

采用ArcGIS对广西数字高程模型数据做预处理，从中提取地理高程数据，在单元节点位置与单元中心位置作曲面差值，得出研究区域范围内的地形数据。对数学模型作地形插值，利用了面积乘积公式：

(11)

图3 局部单元-节点-通道关系Fig.3 Partial unit-node-channel relations

式(11)中：f为等待作插值的点；
fi为已知的点；
Si为以待插值点作顶点的三角形的面积大小；
i∈Sj、i∈Sk分别表示Sj、Sk三角形包含顶点为fi的已知点。

1.5.3 河网处理

图4 河网水系分布Fig.4 Distribution of river network

从基本数据中提取河网信息与水系分布，得到研究区域范围内的河网水系图(图4)。联立控制方程以求解河道水动力信息的方法过于复杂，本文将河网进行数字化处理，基于二叉树河网编码方法进行改进，对河网进行编码分级，为下一步计算做好先后顺序上的准备。

其中，相关定义为：①河道是指没有出流与汇流的河流；
②岔口是指两个及其以上河道的交汇点；
③河段是指组成河道的基本单元；
④过流断面是指河段的衔接处。

编码原则为：将河道岔口的入流端认定为上游，将河道岔口的出流端认定为下游。根据河道、岔口与河段的出现顺序，分别进行编码。若河道的上游岔口只有1条河道，则该河道称为1级编码河道，其上游无河流汇入，下游的河道汇入岔口。每一处岔口，其下游河道级别编码需要等待该岔口上游河道的级别作完编码之后，再继续作编码。按照这个规律往下游作编码，到某个岔口若仅剩唯一河道，则视为该河网水系的出流边界，编码过程如图5所示。根据该编码原则，对研究区域内的双线河网作编码，处理结果共有河道82115段、河道级别最大值已达41，模型局部河网编码如图6所示。

图5 河网编码过程Fig.5 Process of river network coding

2.1 降水因子

降雨对于泥石流的产生，具有明显的推动作用，一方面它是构成泥石流的组分，另一方面它也是大量松散土体运动的载体，进而结合并产生泥石流。本文研究从国家气象局共享数据库中提取了2008年6月8—17日共10日的全国范围内的逐日网格降水量数据，通过自编程序，计算得出模型研究区域内所需的115个节点的降水数据。之所以选择这10日的降雨过程，是因为该降雨过程覆盖了本文全部研究区域且雨量较大，可作为典型暴雨过程。

结合当地典型暴雨分配曲线，将共计10日的日降水过程分配到240 h的时降水过程，以满足数学模型的水动力计算要求。同时，采用泰森多边形法，通过从115个雨量站所采集到的降雨量，推导求出计算范围内的面雨量，从而作进一步的产、汇流与输水过程的数值计算，进而可以得出研究区域范围内网格表面与河网内的积水深度最大值，如图7所示。

2.2 地形因子

模型研究区域范围内的地形条件，毫无疑问是泥石流形成的主要诱发因子之一。本文将坡度作为评价泥石流危险度的指标，利用上文所述的模型地形处理方法，提取等高线、坡度等相关高程信息。之后，根据面积乘积公式，在网格单元节点处与中心处作高程曲面差值，从而得出研究范围内的地形，其三维地形如图8所示，地形高程等值线如图9所示。通过处理得到的地形高程等值线图的走势与广西数字高程模型地形图在大体上保持一致，所以对研究区域范围内地形的实际变化情况有较为真实的反映。

2.3 地质因子

广西当地的地质条件，能够为泥石流的产生提供大量的松散土体。从数据资料中采集整个广西地区的地质分布情况，在数值计算中将研究区域内的砂岩、砾岩、泥岩、页岩等19种不同的地质岩组进

行编码，对这些岩组作数字化处理，从而得到研究范围内各个网格单元上的岩性属性。局部岩组数字化结果如图10所示，将其与来宾市地质分布资料进行比较，说明该地质编码方法具有一定可行性。

图6 河网编码局部放大Fig.6 Partly enlarged and view of river network

图7 模型计算区域网格表面与河网内最大积水深度分布Fig.7 Maximum water depth distribution of mesh surface and river network in model calculation area

图8 模型三维地形Fig.8 3D model terrain

图9 模型地形高程等值线Fig.9 Elevation contours of model terrain

图10 局部地质数字化分布Fig.10 Partial distribution of geological digitization

2.4 植被因子

通过Google Earth软件中的广西卫星地图，获取计算区域内的下垫面植被条件，植被卫星遥感图如图11所示。采用三原色基本原理并结合编程，将红、蓝、绿这三种色素的数值提出，并以绿色数值作为采集植被覆盖率的依据。同时，将绿色数值划分为不同的级别，分别代表林地、草地、田地等，经过数值计算得出研究区域内的植被覆盖数字化信息，如图12所示。比较图11与图12可知，研究区域范围内的植被覆盖在经过具体数据化之后，在大体上能够与卫星遥感数据中的植被覆盖保持相同，表明所采集的植被覆盖率情况是合理可信的。

图11 植被卫星遥感图Fig.11 Remote sensing map of vegetation

图12 植被数字化结果Fig.12 Digitized results of vegetation

2.5 人口情况

将研究区域范围内的人口分布作评价人类活动影响的依据。在数学模型的计算过程中，使用距离加权方法，并与控制面积加权方法作相结合进行使用，把所有行政区域范围内的总人数分拨至各个网格当中，人口分配公式为

(12)

2.6 影响因子权重计算

Pawlak[17]提出了一个关于粗糙集的理论，其特点在于结果完全取决于原始数据，不受主观因素影响。根据粗糙集理论[18-19]，结合基于知识依赖性与信息熵两种条件权重的优点，结合多方面因素分析了条件属性对研究区域中确定、不确定分类子集的影响，采用李大鸣等[20]提出的一种复合的综合客观权重计算方法。使用该方法对各因子权重进行计算，可以客观定量地评定泥石流危险性等级，弥补了先前两种方法的不足，确定各因子权重的方法可以表示为

wi=μwαi+(1-μ)wβi

(13)

式(13)中：wi为客观权值；
wαi为以代数观为依据的权值；
wβi为以信息观为依据的权值；
μ为各个因子。通过计算，得出各个影响因子的权重值，如表1所示。

在模型的数值计算当中，将水流流速、斜坡坡度、地质岩性、人口分布情况与植被覆盖情况分别当作评判泥石流危险度的降雨、地形、地质、人类活动与植被覆盖的指标。在对影响因子作数字化处理的基础上，将各个影响因子的权值乘以其对应的数字化值，从而得到风险值[9]。通过对模型调节，得出一个较为合理的泥石流危险性分级，如表2所示。因此，典型暴雨过程下的泥石流灾害分布，如图13所示。

表1 影响因子权值Table 1 Impact factor parameter

表2 泥石流危险等级Table 2 Debris flow dangerous level

通过国家气象局资源库收集日降雨数据集，选择2011年5月1—10日的降雨数据，计算得出积水深度以及各通道过流量并以其作为水动力推动因素，与斜坡坡度、地质岩性、人口分布情况与植被覆盖情况共同对泥石流危险度做出风险分析。通过泥石流数学模型计算得到该验证时间段内的泥石流灾害危险区域分布，如图14所示；
通过长期统计结果所得到的广西地质灾害分布，如图15所示。

图13 泥石流危险区分布计算结果Fig.13 Calculation results of debris flow hazard zone distribution

图14 验证时间段内泥石流危险区分布计算结果Fig.14 Calculation results of debris flow hazard zone distribution in verification period

图15 地质灾害分布统计结果Fig.15 Statistical results of geological hazard distribution

由图13可得，典型暴雨过程下，广西地区泥石流灾害的高度危险区集中分布于资源县、临桂县、阳朔县、荔浦县、柳江县、柳州市等行政区域，泥石流灾害的危险区集中分布于全州县、龙胜县、灌阳县、融水县、贺州市等行政区域。由图15可得，实际统计结果中，广西地区泥石流多发区集中分布于资源县、临桂县、阳朔县、柳江县等行政区域。因此，通过对比图13与图15可知，实际已发生泥石流灾害的地区，全部位于数学模型计算所得到的泥石流高度危险区范围之内，两者泥石流危险区大体分布基本相同，说明计算结果基本可信。

另外，广西全州县咸水乡洛江村广坑糟屯于2011年5月9日4时20分左右发生了山体滑坡泥石流灾害。由图14可得，在5月1—10日这段验证时间内，整个研究区域只有全州县计算得出了泥石流危险区。并且，受灾地区恰好落在计算所得到的泥石流危险区与高度危险区最为密集的区域范围内(图14右上角红点、蓝块最密集处)，这与计算模拟结果吻合较好，这也进一步验证了计算成果的可信度。

综上所述，无论是长期地质灾害分布统计结果，还是特定时段泥石流灾害事件，都能表明基于水动力学模型的泥石流灾害风险分析具有一定的合理性及可靠性。

以流速、水深作为水动力条件并采用了河网分级编码方法，以坡度、岩性、植被覆盖、人口分布分别作为地形、地质、植被、人口等影响因子的评价指标，采用了一种以粗糙集理论为依据的、更加客观的综合权重计算法，从而构建了泥石流风险分析数学模型。采用广西地区地质灾害分布统计结果与2011年全州县泥石流灾害事件进行模型验证，结果表明，该泥石流风险分析数学模型具有较好的可行性与可信度，在泥石流风险分析的应用当中是合理且可靠的。

针对泥石流破坏力强、规模大、活动频繁等特点，通过水动力模型对泥石流风险进行分析，可以较准确地判断泥石流灾害危险区域分布,为减轻泥石流损害提供指导性技术。这对预防泥石流灾害，减轻泥石流对人类的损害具有重要意义。

猜你喜欢 河网泥石流网格 昆山市平原河网地区活水畅流工程方案设计和效果黑龙江水利科技(2020年8期)2021-01-21追逐作文新天地(初中版)(2019年6期)2019-08-15基于PSR模型的上海地区河网脆弱性探讨华东师范大学学报（自然科学版）(2019年3期)2019-06-24泥石流杂文月刊(2018年21期)2019-01-05基于安卓平台的河网建模与可视化研究电脑知识与技术(2018年7期)2018-06-06重叠网格装配中的一种改进ADT搜索方法北京航空航天大学学报(2017年6期)2017-11-23“民谣泥石流”花粥:唱出自己海峡姐妹(2017年6期)2017-06-24泥石流环球时报(2017-06-14)2017-06-14基于DEM数据的黑河流域信息提取科教导刊·电子版(2016年3期)2016-03-14机械班长科技知识动漫(2016年1期)2016-01-27 相关热词搜索：泥石流，模型，风险，