泥石流堆积扇平面网格化方法与危险性度量方法

摘要

本发明公开了泥石流堆积扇平面网格化方法，以及泥石流堆积扇危险性度量方法。针对现有泥石流堆积扇危险度分区定量评价方法受较多主观性影响的缺陷，本发明首先提供了一种泥石流堆积扇平面网格化方法，该方法在调查数据基础上建立堆积扇平面模型，并对泥石流堆积扇地形划分网格单元。基于该方法，本发明还提供一种泥石流堆积扇危险性度量方法。该方法测算每一网格单元的危险度P(ξ)，并以泥石流堆积扇地形上各网格单元危险度P(ξ)分布特征确定泥石流堆积扇危险性。本发明方法数据来源通过调查、测绘等手段获取，评价结果不受主观因素影响；度量方法将堆积扇危险度分解为P(u)与P(h)两个分量，可以分别度量淤埋与冲击两方面危险性，对泥石流防治工程实际更具有指导性。

说明书

技术领域

本发明涉及一种泥石流堆积扇平面网格化方法，以及以此为基础的 泥石流堆积扇危险性度量方法，属于泥石流灾害防治领域。

背景技术

泥石流堆积扇是泥石流活动的最终产物。在泥石流活动地区，人们 总希望预先知道自己所生活的场所是否安全，以及可能遭遇泥石流泛滥 成灾的危害程度，正如对地震灾害地区进行烈度划分一样，对泥石流活 动区的危险程度等级划分和危险范围确定就成为泥石流灾害地区最重要 的非工程减灾措施之一。

泥石流堆积扇危险度分区是根据影响泥石流堆积泛滥的各因子综合 判定而划分的。目前国内外对泥石流堆积扇的危险性评价方法主要包括 数理分析和数值模拟两大类。前者主要通过建立危险性评价指标体系， 应用数理统计模型进行堆积扇的危险性评价；后者主要通过建立泥石流 的控制方程，应用计算机仿真技术对泥石流堆积扇进行危险性评价。

现有技术《泥石流堆积扇危险度分区定量评价研究》公开了一种泥 石流堆积扇危险度分区度量方法，该方法根据专家综合调查，确定了5个 主导因子参与泥石流危险度评价，并应用德菲尔分析方法计算了各因子 的权重，由此建立了泥石流堆积扇危险度评价模型。该方法是以专家评 分为基础，在建立泥石流堆积扇危险性评价指标体系和确定各个指标权 重的过程中均可能受到主观因素的影响。更重要的是由于人们对于泥石 流的发生机理也并未完全理解，因此，德菲尔分析方法适用过程中的主 观作用更强，可能影响到度量方法的推广性。

发明内容

本发明的目的就是针对现有技术的不足，提供客观实用的泥石流堆 积扇危险性度量方法，该方法在泥石流堆积扇平面模型的网格单元强度 指数测算结果基础上采用网格单元危险度P(ξ)分布特征确定泥石流堆积 扇危险性，方法数据来源通过调查、测量、测绘等手段获取，评价结果 不受主观因素影响。

为实现上述目的，本发明首先提供一种泥石流堆积扇平面网格化方 法，其技术方案如下：

一种泥石流堆积扇平面网格化方法，通过现场调查测绘确定基本地 形数据，利用所述基本地形数据建立泥石流堆积扇平面模型，并对泥石 流堆积扇地形划分网格单元，确定每一网格单元的中心点坐标（x,y）；y 值相同的网格单元确定堆积扇地形上一横断面，根据堆积扇地形特征调 查或测绘获取每一横断面宽度值b；所述泥石流堆积扇平面模型坐标系 以泥石流沟口断面中心为原点o，沿平行沟口断面方向为x轴，沿垂直 沟口断面方向为y轴，沿垂直平面xoy指向原始地面方向为z轴；其特 征在于：所述平面模型依式1建立：

$b=\left(\begin{array}{cc}{b}_0+2y\tan \theta & 0\le y\le \frac{b_m-b_0}{2}\cot \theta \\ \frac{\sqrt{{b_m}^2-{\left(2y\sin \theta {+b}_0\cos \theta \right)}^2}}{\sin \theta }& \frac{b_m-b_0}{2}\cot \theta <y\le l_m\end{array}\right)$   式1

式中，x,y——网格单元的中心点坐标值，

b——网格单元所在的堆积扇横断面的宽度值，m，

b₀——泥石流沟口断面宽度，m，调查测绘确定，

b_m——泥石流堆积扇的最大堆积宽度，m，调查测绘确定 或依经验公式确定，

θ——堆积扇平面所对应的圆周角，°，依式2确定，

$\theta =\arcsin \frac{b_m}{\sqrt{4{l_m}^2+{b_0}^2}}-\arctan \frac{b_o}{2l_m}$0＜θ＜90°   式2

式中，l_m——泥石流堆积扇的最大冲出距离，m，调查测绘确定 或依经验公式计算确定。

上述泥石流堆积扇平面网格化方法的关键在于建立堆积扇平面模 型，利用平面模型确定每一网格单元的中心点坐标（x,y）。其基本技术原 理包括两方面：

第一，平面模型的建立需首先建立地形坐标系。由于泥石流体在堆 积过程中呈直进性，因此建立以泥石流沟口断面中心为原点o，沿平行 沟口断面方向为x轴（横断面），沿垂直沟口断面方向为y轴（纵断面）， 沿垂直平面xoy指向原始地面方向为z轴（z轴的正方向为指向泥石流堆 积前的地面，也即z轴正向朝里）的三维坐标系能够方便建立泥石流堆 积扇的平面模型，以及描述泥石流体在堆积扇上的厚度分布关系和速度 分布关系。xyz坐标系统如图1所示。

第二、泥石流堆积扇的平面模型建立方法的技术原理在于：泥石流 的堆积过程极为复杂，其形成的堆积扇平面常呈现出各种不规则的几何 形态，目前国内外对此均尚未建立起精确统一的堆积模型。经大量野外 调查和几何模型验证发现，实际中，大多数泥石流堆积扇的几何形态可 以近似为扇形平面。由于不同形成概率的堆积扇，其平面范围是不一样 的，其大小主要受泥石流沟口断面宽度b₀、泥石流堆积扇的最大冲出距 离l_m、泥石流堆积扇的最大堆积宽度b_m三个参数影响，因此可以利用此 三个参数来建立泥石流堆积扇的平面模型。在上述地形坐标系中利用扇 形平面几何关系可建立式1表达的泥石流堆积扇的平面模型用以技术原 理的推导应用。

在建立堆积扇平面模型后，即可采用常规方法，如ArcGIS对堆积扇 地形进行栅格化处理，确定每一网格单元的中心点坐标（x,y）。在划分网 格单元后，y值相同的网格单元位于同一个地形横断面上，该横断面的 宽度值b也可相应确定。

以上述泥石流堆积扇平面网格化方法为基础，本发明进一步提供一 种泥石流堆积扇危险性度量方法，其技术方案如下：

一种泥石流堆积扇危险性度量方法，其特征在于：通过现场调查测 绘确定基本地形数据，建立式1表达的堆积扇平面模型；测算每一网格 单元的危险度P(ξ)；以泥石流堆积扇地形上各网格单元危险度P(ξ)分布 特征确定泥石流堆积扇危险性；

所述泥石流堆积扇网格单元危险度P(ξ)包括以速度强度指数表 示的泥石流堆积扇网格单元危险度P(u)与以厚度强度指数表示的泥 石流堆积扇网格单元危险度P(h)，依式3计算确定：

$\left(\begin{array}{c}P\left(u\right)=P_H·I_{\mathrm{uH}}\\ P\left(h\right)=P_H·I_{\mathrm{hH}}\end{array}\right)$   式3

式中，P(u)——以速度强度指数表示的泥石流堆积扇网格单元危险 度，

P(h)——以厚度强度指数表示的泥石流堆积扇网格单元危险 度，

P_H——泥石流堆积扇形成概率，依式4计算确定，

I_uH——速度强度指数，依式5计算确定，

I_hH——厚度强度指数，依式6计算确定，

P_H=1/T   式4

式中，T——降雨周期，年，查询当地水文手册资料确定；

$I_{\mathrm{uH}}=\frac{b-2\left|x\right|}{b}\sqrt{1-\frac{y}{l_m}}$   式5

式中，x,y——网格单元的中心点坐标值，

l_m——泥石流堆积扇的最大冲出距离，m，调查测绘确定或依 经验公式计算确定，

b——网格单元所在的堆积扇横断面的宽度值，m，调查或测 绘确定；

$I_{h_H}=\left(\begin{array}{cc}\exp [-\frac{y^2}{2{C_L}^2{l_m}^2}-\frac{x^2}{2{C_T}^2{(b_0+2y\tan \theta )}^2}]& 0\le y\le \frac{b_m-b_0}{2}\cot \theta \\ \exp [-\frac{y^2}{2{C_L}^2{l_m}^2}-\frac{x^2{\sin }^2\theta }{2{C_T}^2({b_m}^2-{(2y\sin \theta +b_0\cos \theta )}^2)}& \frac{b_m-b_0}{2}\cot \theta <y\le l_m\end{array}\right)$   式6

式中，C_T——堆积扇横断面厚度分布特征常数，依现有技术取值，

C_L——堆积扇纵断面厚度分布特征常数，依现有技术取值；

所述以泥石流堆积扇地形上各网格单元危险度P(ξ)分布特征确定泥 石流堆积扇危险性的判别条件如表1所示。

表1泥石流堆积扇危险性判别条件

上述泥石流堆积扇危险性度量方法，是在测算泥石流堆积扇上划分 的每一网格单元在堆积扇形成概率P_H时的危险度P(ξ)的基础上，结合扇 形成概率P_H与危险度P(ξ)判断泥石流堆积扇不同区域的危险性，其技术 思想与原理在于：

第一，泥石流堆积扇的危险度主要包括两层含义：一是给定规模堆 积扇的形成概率（假定堆积扇的形成概率就是对应泥石流的发生概率）； 二是给定规模堆积扇上各单元区内所表现的破坏强度（即网格单元的强 度指数），因此泥石流堆积扇的危险度P(ξ)可以由堆积扇的形成概率及堆 积扇的强度指数两部分表达，即式7：

$P\left(\xi \right)=P_H·I_{{\xi }_H}$   式7

上式中，表示形成概率为P_H的堆积扇上变量ξ描述的强度指数。

第二，关于泥石流堆积扇形成概率P_H。泥石流堆积扇形成概率的确 定，采用泥石流发生概率的确定方法。泥石流发生概率通常采用年超越 概率（即t＝1年内发生次数≥1的概率）表示，假定单个泥石流为独立的 随机点事件，则其发生的年超越概率可直接运用研究自然界中随机点事 件的概率模型进行计算，一般地，用于研究自然界中的随机点事件概率 模型有Poisson模型和Binomial模型，前者是连续时间模型，后者是离 散时间模型。本发明方法中采用Binomial模型计算泥石流发生概率，当 泥石流的发生概率无法直接确定时，假定泥石流与降雨同频率、且同步 发生，则有P_H=1/T。式中T表示降雨周期，均可由查询当地水文手册资 料获取。

第三，关于泥石流堆积扇强度指数给定概率的堆积扇破坏强度 可以用堆积扇上各网格单元的变量ξ来描述，其表达式如式8：

$I_{{\xi }_H}=\frac{\xi (x,y)}{{\xi }_m}$   式8

上式中，ξ(x,y)——堆积扇平面上任一点的ξ值，ξ_m——堆积扇平面上ξ的 最大值（假定堆积扇平面内的最大堆积厚度与最大流体速度均位于坐标 原点）。对于堆积扇的危险性测算通常需要考虑流体堆积厚度h与速度u 两个变量，由此对应强度指数的两个分量，即速度强度指数I_uH与厚度 强度指数I_hH。当确定了堆积扇上的最大堆积厚度与最大流体速度值时， 根据强度指数的定义式8，由网格单元的厚度（或速度）强度指数值可 以求出该网格单元上的厚度（或速度）值。

进而，对于堆积扇上的每一个网格单元的危险度P(ξ)都具有两个分 量，即以速度强度指数表示的泥石流堆积扇网格单元危险度P(ξ＝u)（简 写作P(u)）与以厚度强度指数表示的泥石流堆积扇网格单元危险度 P(ξ＝h)（简写作P(h)）。在度量泥石流堆积扇危险性时，这两个分量都 应该被测算。

第四，关于网格单元厚度强度指数I_hH。根据Yuan-Fan Tsai的实验结 论，泥石流堆积扇在其横断面y＝y_c上的厚度分布近似满足高斯分布曲线 其中z—堆积扇上点(x,y_c)的厚度值(m)，z_c—横断 面y＝y_c在y轴上的厚度值(m)；b—横断面y＝y_c的宽度值(m)；C_T—堆积 扇横断面厚度分布特征常数。泥石流堆积扇在其纵断面x＝0上的厚度分 布近似满足半高斯分布曲线其中z_c—堆积扇上点 (0,y)的厚度值(m)，z_m—y轴上的最大厚度值（即堆积扇上的最大厚度 值h_m，处于出山口位置）(m)，l_m—堆积扇的最大冲出距离(m)，C_L—堆 积扇纵断面厚度分布特征常数。此二分布曲线方程用于表示每个网格单 元中心点的Z轴数值，即堆积扇在该点的淤埋厚度分布情况。联立此二 分布曲线方程与堆积扇平面关系式（即式1），可以推导得堆积扇上各网 格单元的厚度强度指数I_hH表达式，即式6。

根据Yuan-Fan Tsai的实验结论，堆积扇横断面厚度分布特征常数C_T取值范围0.204～0.225，该值大小主要取决于泥石流体的颗粒级配特征和 流体的稠度。具体可根据综合判定泥石流的性质进行取值，当泥石流近 似为水石型时，两者均取最大值，随着泥石流稠度增大，取值降低。堆 积扇纵断面厚度分布特征常数C_L取值范围0.373～0.401，其具体取值方法 与C_T值相同。参见参考文献1。

特殊地，当θ＝90°时，式6则变形为式9：

$I_{h_H}=\exp [-\frac{2y^2}{{C_L}^2{b_m}^2}-\frac{x^2}{2{C_T}^2({b_m}^2-4y^2)}]$   式9

第五，关于网格单元速度强度指数I_uH。影响泥石流堆积扇上流速衰 减的因素众多，本方法中假设这些影响因素对堆积扇主流线上流速衰减 所起的力学效果在空间上是恒定的，则泥石流在主流线上的堆积过程可 以简化为匀减速运动；同时，假定流速在横断面上的衰减过程与距主流 线的距离呈线性关系，且有堆积扇边界上的速度值为零，则可以推导得 堆积扇上各网格单元的速度强度指数I_uH表达式，即式5。

第六，结合堆积扇形成概率P_H与危险度P(ξ)作为泥石流堆积扇各网 格单元危险等级的判别条件，其技术原理在于：当堆积扇的形成概率P_H越 小时，堆积扇的规模就越大，其破坏程度相应愈强，即在P_H较小的情况 下，P(ξ)值很小的网格单元的危险等级也可能较高，因此，在进行单元 危险等级划分时，需要综合考虑P_H和P(ξ)两者因素。最终根据不同危险 等级所占的堆积扇面积的比例范围可以度量堆积扇的淤埋或冲击能力。

上述泥石流堆积扇危险性度量方法中，对堆积扇地形栅格化处理时， 网格单元的大小应当依据评价精度要求确定。例如，若地形图上房屋的 平面为20cm×10cm，则单元网格就应尽量小于20cm×20cm。

本方法适用时，对于已经形成过的堆积扇，可以直接利用野外测量 测绘获取基础数据，但对于未形成过或没有历史数据的堆积扇l_m与b_m无 法通过野外测量测绘获取，则需要选取合适的用于估算泥石流堆积扇最 大冲出距离和最大堆积宽度的经验公式计算确定。如汶川震区泥石流堆 积扇的最大冲出距离计算式：l_m＝0.36A^0.06+0.03(V_L·H)^0.54-0.18，最大堆积 宽度计算式：b_m＝0.40A^0.08+0.04(V_L·H)^0.35-0.23，其中，A为泥石流沟的流 域面积（km²），V_L为流域的动储量（10⁶m³），H为流域的高差（km）， 三者均可以通过遥感影像解译或野外调查获取。

本发明方法应用于泥石流的堆积区有形成扇形形态的地形条件。对 于不满足该条件的泥石流堆积扇危险性度量，可以首先根据当地地形条 件适当修正本方法中提出的平面模型，再基于本度量方法的技术原理进 行危险性度量。

与现有技术相比，本发明的有益效果是：（1）本发明方法在泥石流 堆积扇平面模型的网格单元强度指数测算结果基础上采用网格单元危险 度P(ξ)分布特征确定泥石流堆积扇危险性，方法数据来源通过调查、测 量、测绘等手段获取，评价结果不受主观因素影响。；（2）本发明方法将 堆积扇每一位点的危险度分解为P(u)与P(h)两个分量，前者用于评价堆 积扇的冲击破坏危险性，后者用于评价堆积扇的淤埋破坏危险性，对堆 积扇的危险度的度量更为具体，对泥石流防治工程实践有更具体明确的 解读价值。

参考文献：Yuan-Fan Tsai.Three-Dimensional Topography of  Debris-Flow Fan.Journal of Hydraulic Engineering(2006),132:307-318.

附图说明

图1是堆积扇的坐标系统图。

图2是技术流程示意图。

图3麻柳沟泥石流堆积扇基于堆积厚度的危险分区图。

图4麻柳沟泥石流堆积扇基于流体速度的危险分区图。

具体实施方式

下面结合附图，对本发明的优选实施例作进一步的描述。

实施例一

如图1～图4所示，用本发明方法对四川省都江堰市龙池镇麻柳沟 泥石流堆积扇危险性进行度量。图2是本实施例技术流程示意图。

2010年8月14日，四川省都江堰市龙池镇麻柳沟暴发泥石流，在 其堆积扇范围内淤埋农房3户、集中安置区安置房1栋。堆积扇上流体 速度最大值为6.48m/s，堆积厚度最大值为5.5m。

1、实地调查测绘

经实地调查测绘，确定基本参数数据：麻柳沟泥石流沟口断面宽度 b₀=30m，泥石流堆积扇的最大冲出距离l_m=218m，泥石流堆积扇的最 大堆积宽度b_m=210m。

如图1所示，以泥石流沟口断面中心为原点o，沿平行沟口断面方 向为x轴，沿垂直沟口断面方向为y轴，沿垂直xoy平面指向原始地面 为z轴。

2、确定泥石流堆积扇形成概率P_H

查询四川省水文手册资料确定触发该次泥石流的降雨为50年一遇， 即降雨周期T＝50，代入式2则有堆积扇形成概率P_H=1/50=0.02。

3、建立堆积扇平面模型

将参数b₀、l_m、b_m代入式2计算确定泥石流堆积扇平面所对应的圆 周角θ，得到θ＝24.8°。

将参数b₀、l_m、b_m、θ代入式1建立堆积扇平面模型，得到下式：

$b=\left(\begin{array}{cc}30+0.92y& 0\le y\le 194.93\\ \frac{\sqrt{44100-{(0.84y+27.24)}^2}}{0.42}& 194.93<y\le 218\end{array}\right)$

4、地形栅格化处理

根据堆积扇平面模型对堆积扇平面进行网格单元划分，对麻柳沟泥 石流堆积扇平面进行5m×5m的网格单元划分，整个堆积区共划分出1152 个网格单元，确定每一网格单元的中心点坐标（x,y）。

5、计算网格单元速度强度指数I_uH、厚度强度指数I_hH

将参数l_m、b及各中心点坐标（x,y）代入式5，计算确定每一网格单 元中心点的速度强度指数I_uH，部分数据如表2。

表2部分网格单元的计算结果值

序号 x/m y/m b/m I_uHI_hHP(u) P(h) 1 -1.22 0.21 30 0.9189 0.9999 0.0184 0.0200 2 -4.20 5.21 35 0.7494 0.9972 0.0150 0.0199 3 -7.08 11.67 41 0.6348 0.9877 0.0127 0.0198 4 19.48 12.25 41 0.0543 0.9722 0.0011 0.0194 5 12.30 20.30 49 0.4710 0.9631 0.0094 0.0193 6 13.00 36.10 63 0.5376 0.9003 0.0108 0.0180 7 -36.00 45.24 72 0.0046 0.8067 0.0001 0.0161 8 -22.55 55.09 81 0.3813 0.7778 0.0076 0.0156 9 19.08 65.90 91 0.4835 0.7102 0.0097 0.0142 10 -35.85 76.26 100 0.2291 0.6051 0.0046 0.0121 11 -11.80 83.10 106 0.6122 0.5933 0.0122 0.0119 12 56.75 94.29 117 0.0209 0.4267 0.0004 0.0085 13 67.94 102.16 124 0.0699 0.3458 0.0014 0.0069 14 71.55 110.23 131 0.0626 0.2898 0.0013 0.0058 15 -53.11 115.02 136 0.1497 0.3086 0.0030 0.0062 16 45.24 120.13 141 0.2386 0.2957 0.0048 0.0059 17 47.74 140.59 159 0.2388 0.1888 0.0048 0.0038 18 8.60 150.25 168 0.5005 0.1848 0.0100 0.0037 19 88.28 162.52 180 0.0083 0.0566 0.0002 0.0011 20 -92.31 194.32 209 0.0381 0.0105 0.0008 0.0002

根据野外调查确定麻柳沟泥石流为粘性泥石流，测试确定容重γ约 为1.92t/m³，故常数C_T取0.220，C_L取0.390，将参数l_m、C_T、C_L、b₀、θ 及各中心点坐标（x,y）代入式6，计算确定每一中心点的厚度强度指数I_hH， 部分数据如表2。

6、计算网格单元危险度P(ξ)

将每一网格单元的参数P_H、I_uH、I_hH分别代入式3，计算每一网格单 元的P(u)与P(h)，部分数据如表2。

7、判断泥石流堆积扇危险性

根据表1泥石流堆积扇危险性判别条件，判断堆积扇各网格单元的 危险性，P(u)用于评价堆积扇的冲击破坏危险性，P(h)用于评价堆积扇 的淤埋破坏危险性。分别绘制麻柳沟泥石流堆积扇平面范围的危险分区 图，如图3、图4所示。

图3是基于堆积厚度的危险分区图。图3中，P(h)极高危险性区域 高危险性区域、中危险性区域、低危险性区域所占的比例分别为26.2%， 19.0%，36.0%，18.8%。分区结果表明，泥石流形成后，堆积扇的厚度 强度指数大于0.25（根据堆积厚度最大值为5.5m，即堆积扇厚度大于0.25 ×5.5=1.375m）的范围占45.2%，表现出很强的淤埋能力，与实际情况 大致吻合。

图4是基于流体速度的危险分区图。图4中，P(u)极高危险性区域 高危险性区域、中危险性区域、低危险性区域所占的比例分别为16.7%， 31.3%，37.8%，14.2%。分区结果表明，泥石流形成后，堆积扇的速度 强度指数居于0.05～0.5之间（利用堆积扇上流体速度最大值6.48m/s， 计算得到即泥石流速度在0.05×6.48m＝0.324至0.5×6.48m＝3.24m/s） 的范围占69.1%，表现出较强的冲击能力。