本发明涉及一种泥石流坡面物源启动量动态计算方法及系统，属于泥石流灾害防控领域，将待计算泥石流物源坡面划分为土柱后，确定选定的中心土柱的周围6个相邻土柱的位置关系，然后利用极限分析上限定理计算土柱的最不利滑动面及滑动面上的不平衡力，进而判断最不利滑动面是否发生失稳，根据中心土柱的侧向拉应力的连接键断裂情况，计算失稳土柱与其周围土柱的作用力方式和大小，最终判断土柱是否流态化，得到泥石流坡面物源启动量。本发明能对潜在的坡面泥石流的启动流量进行准确计算，将有利于解决长期存在的泥石流启动量计算难、精确度低的问题。
1.一种泥石流坡面物源启动量动态计算方法，其特征在于，所述计算方法包括：选择待计算的泥石流物源坡面，将所述待计算的泥石流物源坡面划分为规则排列的i行j列的正六边形的土柱；
根据极限分析上限定理计算每个土柱的最不利滑动面，并计算最不利滑动面上的不平衡力；
所述根据极限分析上限定理计算每个土柱的最不利滑动面，并计算最不利滑动面上的不平衡力，具体包括：
根据极限分析上限定理，计算最不利滑动面的深度Hsd和角度α+β：总的外功率为：
计算最不利滑动面上的不平衡力，具体包括：根据 计算
利用所述最不利滑动面的下滑力和所述最不利滑动面的抗滑力计算所述最不利滑动面上的不平衡力；
其中： τT为周围相邻6个土柱对中心土柱的应力总和，τ′Tn为所述周围相邻土柱对所述中心土柱的应力；
Lg表示土柱间的间距，c为待计算的泥石流物源根土复合体的粘聚力，Hsd为最不利滑动面的深度；θ为所述待计算泥石流物源坡面土体含水率，ρw为水的密度， 为所述待计算泥石流物源坡面根土复合体的孔隙率，ρr为所述待计算泥石流物源坡面根土复合体的密度，g为重力加速度，α为基岩的角度，β为最不利滑动面同基岩的角度，χ为与水土特征曲线有关的参数，h为毛细管压力水头，γ为待计算的泥石流物源坡面根土复合体的内摩擦角，τrf为待计算泥石流物源坡面径流侵蚀力，γw为水的重度，hs为待计算泥石流物源坡面所在区域的潜水位水头，J为待计算泥石流物源坡面所在区域的潜水位水力坡度，H1为所述中心土柱的坡面高程，H2为所述中心土柱的基岩高程，Hd为待计算泥石流物源坡面所在区域的地下水位；
根据所述最不利滑动面上的不平衡力分别判断相应土柱是否发生失稳；若判断结果为否，则该土柱停止计算；若判断结果为是，则根据土柱是否失稳，以及中心土柱与其周围六个相邻土柱的重心高程，确定所述中心土柱失稳方向，并获取周围六个相邻土柱对中心土柱的作用力方式；
利用纤维束模型模拟土柱间的相互作用，并根据连接键断裂数量计算中心土柱与相邻土柱间的应力大小；若判断结果为侧向拉应力的连接键未全部断裂，则该土柱停止计算；若判断结果为侧向拉应力的连接键全部断裂，则根据流态化准则判断土柱是否流态化，若判断结果为否，则该土柱停止计算；若判断结果为是，则该土柱启动，记录坡面物源泥石流启动量。
2.根据权利要求1所述的泥石流坡面物源启动量动态计算方法，其特征在于，所述待计算的泥石流物源根土复合体的粘聚力c与所述待计算的泥石流物源坡面根土复合体的内摩擦角γ服从威布尔分布：
所述待计算泥石流物源坡面土体含水率θ、所述待计算泥石流物源坡面根土复合体的孔隙率 与待计算泥石流物源坡面根土复合体的密度ρ服从正态分布：
式中：μθ, μρ分别表示含水率、孔隙率和密度服从的正态分布的期望，σθ, σρ分别表示含水率、孔隙率和密度服从的正态分布的方差，期望和方差根据坡面取得的实测值确定。
3.根据权利要求1所述的泥石流坡面物源启动量动态计算方法，其特征在于，所述根据土柱是否失稳，以及中心土柱与其周围六个相邻土柱的重心高程，确定所述中心土柱失稳方向，并获取周围六个相邻土柱对中心土柱的作用力方式，具体包括：获取所述中心土柱的位置(i,j)，并确定所述中心土柱的重心坐标为(xi，j,yi,j,zi,j)；
根据位置坐标和几何原理得到所述周围相邻土柱的重心坐标(xk,yk,zk)(k＝0,1,2,3,
根据所述中心土柱的重心坐标和周围相邻土柱的重心坐标得到周围相邻土柱相对中心土柱的向量(x'k,y'k,z'k)， Lg为相邻土柱间的距离；
周围相邻土柱中，两相邻土柱的和向量为(x′k,k+1,y′k,k+1,z′k,k+1)所述和向量投影在二维平面内的向量为Dk,k+1(k＝0,1,2,3,4,5)；
所述中心土柱运动方向在二维平面投影向量(x″i,j,y″i,j)；
根据最小作用量原理，取重心最低min(z′k,k+1)(k＝0,1,2,3,4,5)的和向量为土柱运动的方向向量；
k取min(z′k,k+1)(k＝0,1,2,3,4,5)时的值，所述中心土柱在运动方向在二维平面的投影向量为(x″i,j,y″i,j)：x″i,j＝x′k,k+1、y″i,j＝y′k,k+1；
周围土柱运动方向在二维平面投影向量(x″k,y″k)；
当相邻土柱的运动方向向量的夹角为锐角时，土柱间的作用力为压应力，当相邻土柱间的运动方向向量为钝角的时候，土柱间的作用力为拉应力，若 则：τ′Tk＝‑τTk；
τTk为周围土柱失稳时对中心土柱的作用力大小，τ′Tk为中心土柱计算下滑力时受到的周围土柱的力。
4.根据权利要求1所述的泥石流坡面物源启动量动态计算方法，其特征在于，所述利用纤维束模型模拟土柱间的相互作用，并根据连接键断裂数量计算中心土柱与相邻土柱间的应力大小，具体包括：
采用纤维束模型模拟土柱间的相互作用，根据纤维束的破坏情况以及土柱间的接触面积，计算土柱破坏过程应力分配的大小；
同理可得，τT1，τT2，τT3，τT5的计算公式，式中n＝0,1,2,3,4,5；
式中：Acn表示的是土柱与周围六个的土柱间的接触面积，τT0、τT1、τT2、τT3表示的是i,j土柱对周围土柱的拉应力，τT4、τT5表示的是i,j土柱对周围土柱的压应力，Nn表示i,j土柱与编号为n的土柱间的完整纤维的数量，K2/K1表示压应力与拉应力间的分配比，hrf表示中心土柱的侵蚀深度，hrfn表示周围土柱编号为n的土柱的侵蚀深度，H1表示中心土柱坡面的高程，AH为土柱的截面面积，Lg为相邻土柱间的距离；H1n表示周围土柱编号为n的土柱的坡面高程，τR为土柱沿最不利滑动面破坏后的残余强度，Hsd为最不利滑动面的深度。
5.一种泥石流坡面物源启动量动态计算系统，其特征在于，所述计算系统包括：坡面划分模块，用于选择待计算的泥石流物源坡面，将所述待计算的泥石流物源坡面划分为规则排列的i行j列的正六边形的土柱；
滑动面及不平衡力计算模块，用于根据极限分析上限定理计算每个土柱的最不利滑动面，并计算最不利滑动面上的不平衡力；
所述滑动面及不平衡力计算模块具体包括：最不利滑动面深度和角度计算单元，用于根据极限分析上限定理，计算最不利滑动面的深度Hsd和角度α+β：
最不利滑动面上的不平衡力计算单元，用于计算最不利滑动面上的不平衡力，具体包括：
最不利滑动面上的不平衡力计算子单元，用于利用所述最不利滑动面的下滑力和所述最不利滑动面的抗滑力计算所述最不利滑动面上的不平衡力；
其中： τT为周围相邻6个土柱对中心土柱的应力总和，τ′Tn为所述周围相邻土柱对所述中心土柱的应力；
Lg表示土柱间的间距，c为待计算的泥石流物源根土复合体的粘聚力，Hsd为最不利滑动面的深度；θ为所述待计算泥石流物源坡面土体含水率，ρw为水的密度， 为所述待计算泥石流物源坡面根土复合体的孔隙率，ρr为所述待计算泥石流物源坡面根土复合体的密度，g为重力加速度，α为基岩的角度，β为最不利滑动面同基岩的角度，χ为与水土特征曲线有关的参数，h为毛细管压力水头，γ为待计算的泥石流物源坡面根土复合体的内摩擦角，τrf为待计算泥石流物源坡面径流侵蚀力，γw为水的重度，hs为待计算泥石流物源坡面所在区域的潜水位水头，J为待计算泥石流物源坡面所在区域的潜水位水力坡度，H1为所述中心土柱的坡面高程，H2为所述中心土柱的基岩高程，Hd为待计算泥石流物源坡面所在区域的地下水位；
失稳判断模块，用于根据所述最不利滑动面上的不平衡力分别判断相应土柱是否发生失稳；若判断结果为否，则该土柱停止计算；若判断结果为是，则失稳方向判定模块，用于根据土柱是否失稳，以及中心土柱与其周围六个相邻土柱的重心高程，确定所述中心土柱失稳方向，并获取周围六个相邻土柱对中心土柱的作用力方式；
连接键断裂判断模块，用于利用纤维束模型模拟土柱间的相互作用，并根据连接键断裂数量计算中心土柱与相邻土柱间的应力大小；若判断结果为侧向拉应力的连接键未全部断裂，则该土柱停止计算；若判断结果为侧向拉应力的连接键全部断裂，则流态化判断模块，用于根据流态化准则判断土柱是否流态化，若判断结果为否，则该土柱停止计算；若判断结果为是，则记录模块，用于记录该土柱启动，记录坡面物源泥石流启动量。
6.根据权利要求5所述的泥石流坡面物源启动量动态计算系统，其特征在于，所述待计算的泥石流物源根土复合体的粘聚力c与所述待计算的泥石流物源坡面根土复合体的内摩擦角γ服从威布尔分布：
所述待计算泥石流物源坡面土体含水率θ、所述待计算泥石流物源坡面根土复合体的孔隙率 与待计算泥石流物源坡面根土复合体的密度ρ服从正态分布：
式中：μθ, μρ分别表示含水率、孔隙率和密度服从的正态分布的期望，σθ, σρ分别表示含水率、孔隙率和密度服从的正态分布的方差，期望和方差根据坡面取得的实测值确定。
7.根据权利要求5所述的泥石流坡面物源启动量动态计算系统，其特征在于，所述失稳方向判定模块具体包括：
中心土柱位置获取单元，用于获取所述中心土柱的位置(i,j)，并确定所述中心土柱的重心坐标为(xi，j,yi,j,zi,j)；
周围相邻土柱的重心坐标计算单元，用于根据位置坐标和几何原理得到所述周围相邻土柱的重心坐标(xk,yk,zk)(k＝0,1,2,3,4,5,6)为作用向量计算单元，用于根据所述中心土柱的重心坐标和周围相邻土柱的重心坐标得到周围相邻土柱相对中心土柱的向量(x'k,y'k,z'k)，Lg为相邻土柱间的距离；
周围相邻土柱中，两相邻土柱的和向量为(x′k,k+1,y′k,k+1,z′k,k+1)所述和向量投影在二维平面内的向量为Dk,k+1(k＝0,1,2,3,4,5)；
所述中心土柱运动方向在二维平面投影向量(x″i,j,y″i,j)；
土柱运动的方向向量确定单元，用于根据最小作用量原理，取重心最低min(z′k,k+1)(k＝0,1,2,3,4,5)的和向量为土柱运动的方向向量；
k取min(z′k,k+1)(k＝0,1,2,3,4,5)时的值，所述中心土柱在运动方向在二维平面的投影向量为(x″i,j,y″i,j)：x″i,j＝x′k,k+1、y″i,j＝y′k,k+1；
周围土柱运动方向在二维平面投影向量(x″k,y″k)；
当相邻土柱的运动方向向量的夹角为锐角时，土柱间的作用力为压应力，当相邻土柱间的运动方向向量为钝角的时候，土柱间的作用力为拉应力，若 则：τ′Tk＝‑τTk；
τTk为周围土柱失稳时对中心土柱的作用力大小，τ′Tk为中心土柱计算下滑力时受到的周围土柱的力。
8.根据权利要求5所述的泥石流坡面物源启动量动态计算系统，其特征在于，采用纤维束模型模拟土柱间的相互作用，根据纤维束的破坏情况以及土柱间的接触面积，计算土柱破坏过程应力分配的大小；
同理可得，τT1，τT2，τT3，τT5的计算公式，式中n＝0,1,2,3,4,5；
式中：Acn表示的是土柱与周围六个的土柱间的接触面积，τT0、τT1、τT2、τT3表示的是i,j土柱对周围土柱的拉应力，τT4、τT5表示的是i,j土柱对周围土柱的压应力，Nn表示i,j土柱与编号为n的土柱间的完整纤维的数量，K2/K1表示压应力与拉应力间的分配比，hrf表示中心土柱的侵蚀深度，hrfn表示周围土柱编号为n的土柱的侵蚀深度，H1表示中心土柱坡面的高程，AH为土柱的截面面积，Lg为相邻土柱间的距离；H1n表示周围土柱编号为n的土柱的坡面高程，τR为土柱沿最不利滑动面破坏后的残余强度，Hsd为最不利滑动面的深度。
一种泥石流坡面物源启动量动态计算方法及系统
技术领域
[0001] 本发明涉及泥石流灾害防治技术领域，特别是涉及一种泥石流坡面物源启动量动态计算方法与系统。
背景技术
[0002] 近些年，随着全球气候变暖，受强降雨、暴风等极端天气和大地震活跃等因素的影响，泥石流灾害频发，严重影响山区人民的生命财产安全，并制约着我国山区国民经济发展
和社会的可持续发展。泥石流坡面物源启动量的准确计算直接影响泥石流预警预报的准确
性以及防治工程设计的可靠性，其难点在于如何利用降雨数据，并结合泥石流坡面物源启
动物理力学机制对泥石流坡面物源启动量随时间变化进行连续计算。泥石流坡面物源在降
雨作用下，启动过程通常受渗流作用、侵蚀作用、渗流和侵蚀的共同作用的影响。由于影响
因素复杂多样，泥石流具有突发性、随机性、难预测性等特点，造成泥石流坡面物源启动量
难以准确计算。目前，国内外关于泥石流坡面物源启动量的计算方法，多数都局限于经验模
型和基于大量简化的方法，与实际泥石流坡面物源启动量及其启动过程特征差异较大，导
致泥石流启动量估算误差大。
[0003] 现有泥石流坡面物源启动量计算方法存在的问题或缺陷。
[0004] 泥石流坡面物源启动量计算主要基于以下四种模型：随机过程模型、无限边坡模型、坡面侵蚀模型和土柱模型，这四种模型分别反映了泥石流坡面物源的随机启动、重力和
渗流作用、侵蚀作用以及渐进性破坏，但目前泥石流坡面物源启动量计算过程都较为单一。
基于单一模型的计算主要存在的问题如下：
[0005] 1、基于随机过程模型的计算，考虑了降雨、径流冲刷、降雨入渗、重力作用和土体参数的随机分布，但主要是基于统计公式，缺乏物理意义和确定性，计算结果波动较大。
[0006] 2、基于无限边坡模型的计算，将坡面划分成若干条块，基于极限平衡计算条块下滑力和抗滑力。但认为泥石流坡面物源一次性启动，且启动深度与物源厚度一致，这与实际
泥石流坡面物源渐进启动过程不符，也没有考虑径流冲刷、土体参数随机分布和降雨入渗。
[0007] 3、基于坡面侵蚀模型的计算，主要考虑降雨导致的坡面径流对泥石流坡面物源的侵蚀作用，没有考虑降雨入渗、重力作用和土体参数的随机分布。
[0008] 4、基于土柱模型的计算，将坡面离散为土柱单元，考虑了由于重力引起的土柱间的相互作用，但认为每个土柱代表的泥石流坡面物源一次性启动，且坡面启动深度与物源
厚度一致，没有考虑径流冲刷和土体参数随机分布等因素导致的坡面物源渐进启动过程。
[0009] 除了基于以上4种模型的坡面物源的计算，还有一些经验公式或简化方法，但都受到了各种条件的限制，导致泥石流坡面物源启动量计算结果与实际相差太大，难以满足泥
石流预警及防治中对泥石流坡面物源启动量的计算要求。
发明内容
[0010] 本发明的目的是提供一种泥石流坡面物源启动量动态计算方法与系统，可以准确计算泥石流的启动量。
[0011] 为实现上述目的，本发明提供了如下方案：
[0012] 一种泥石流坡面物源启动量动态计算方法，所述计算方法包括：
[0013] 选择待计算的泥石流物源坡面，将所述待计算的泥石流物源坡面划分为规则排列的i行j列的正六边形的土柱；
[0014] 根据极限分析上限定理计算每个土柱的最不利滑动面，并计算最不利滑动面上的不平衡力；
[0015] 根据所述最不利滑动面上的不平衡力分别判断相应土柱是否发生失稳；若判断结果为否，则该土柱停止计算；若判断结果为是，则
[0016] 根据土柱是否失稳，以及中心土柱与其周围六个相邻土柱的重心高程，确定所述中心土柱失稳方向，并获取周围六个相邻土柱对中心土柱的作用力方式；
[0017] 利用纤维束模型模拟土柱间的相互作用，并根据连接键断裂数量计算中心土柱与相邻土柱间的应力大小；若判断结果为侧向拉应力的连接键未全部断裂，则该土柱停止计
算；若判断结果为侧向拉应力的连接键全部断裂，则
[0018] 根据流态化准则判断土柱是否流态化，若判断结果为否，则该土柱停止计算；若判断结果为是，则该土柱启动，记录坡面物源泥石流启动量。
[0019] 可选的，所述根据极限分析上限定理计算每个土柱的最不利滑动面，并计算最不利滑动面上的不平衡力，具体包括：
[0020] 根据极限分析上限定理，计算最不利滑动面的深度Hsd和角度α+β：
[0021] 总的外功率为：
[0022] 滑动面上的内能耗散为：
[0023] 土柱安全系数K可以表示为：滑动面对应着最小的安全系数：
[0024] 计算最不利滑动面上的不平衡力，具体包括：
[0025] 根据 计算最不利滑动面的下滑力；
[0026] 根据计算最不利滑动面的抗滑力；
[0027] 利用所述最不利滑动面的下滑力和所述最不利滑动面的抗滑力计算所述最不利滑动面上的不平衡力；
[0028] 其中： τT为周围相邻6个土柱对中心土柱的应力总和，τ′Tn为所述周围相邻土柱对所述中心土柱的应力；
[0029] Lg表示土柱间的间距，c为待计算的泥石流物源根土复合体的粘聚力，Hsd为最不利滑动面的深度；θ为所述待计算泥石流物源坡面土体含水率，ρw为水的密度， 为所述待计
算泥石流物源坡面根土复合体的孔隙率，ρr为所述待计算泥石流物源坡面根土复合体的密
度，g为重力加速度，α为基岩的角度，β为最不利滑动面同基岩的角度，χ为与水土特征曲线
有关的参数，h为毛细管压力水头，γ为待计算的泥石流物源坡面根土复合体的内摩擦角，
τrf为待计算泥石流物源坡面径流侵蚀力，γw为水的重度，hs为待计算泥石流物源坡面所在
区域的潜水位水头，J为待计算泥石流物源坡面所在区域的潜水位水力坡度，H1为所述中心
土柱的坡面高程，H2为所述中心土柱的基岩高程，Hd为待计算泥石流物源坡面所在区域的地
下水位。
[0030] 可选的，所述待计算的泥石流物源根土复合体的粘聚力c与所述待计算的泥石流物源坡面根土复合体的内摩擦角γ服从威布尔分布：
[0031]
[0032] 所述待计算泥石流物源坡面土体含水率θ、所述待计算泥石流物源坡面根土复合体的孔隙率 与待计算泥石流物源坡面根土复合体的密度ρ服从正态分布：
[0033]
[0034] 式中：μθ, μρ分别表示含水率、孔隙率和密度服从的正态分布的期望，σθ2,2
σρ分别表示含水率、孔隙率和密度服从的正态分布的方差，期望和方差根据坡面取得的实
测值确定。
[0035] 可选的，所述根据土柱是否失稳，以及中心土柱与其周围六个相邻土柱的重心高程，确定所述中心土柱失稳方向，并获取周围六个相邻土柱对中心土柱的作用力方式，具体
包括：
[0036] 获取所述中心土柱的位置(i,j)，并确定所述中心土柱的重心坐标为(xi，j,yi,j,zi,j)；
[0037] 根据位置坐标和几何原理得到所述周围相邻土柱的重心坐标(xk,yk,zk)(k＝0,1,2,3,4,5,6)为
[0038]
[0039] 根据所述中心土柱的重心坐标和周围相邻土柱的重心坐标得到周围相邻土柱相对中心土柱的向量(x'k,y'k,z'k)， Lg为相邻土柱间
的距离；
[0040] 周围相邻土柱中，两相邻土柱的和向量为(x′k,k+1,y′k,k+1,z′k,k+1)
[0041]
[0042] 所述和向量投影在二维平面内的向量为Dk,k+1(k＝0,1,2,3,4,5)；
[0043] 所述中心土柱运动方向在二维平面投影向量(x″i,j,y″i,j)；
[0044] 根据最小作用量原理，取重心最低min(z′k,k+1)(k＝0,1,2,3,4,5)的和向量为土柱运动的方向向量；
[0045] k取min(z′k,k+1)(k＝0,1,2,3,4,5)时的值，所述中心土柱在运动方向在二维平面的投影向量为(x″i,j,y″i,j)：x″i,j＝x′k,k+1、y″i,j＝y′k,k+1；
[0046] 周围土柱运动方向在二维平面投影向量(x″k,y″k)；
[0047] 相邻土柱间的夹角θk为
[0048] 当相邻土柱的运动方向向量的夹角为锐角时，土柱间的作用力为压应力，当相邻土柱间的运动方向向量为钝角的时候，土柱间的作用力为拉应力，
[0049] 若 则：τ′Tk＝‑τTk；
[0050] 若 则：τ′Tk＝0；
[0051] 若 则：τ′Tk＝τTk；
[0052] τTk为周围土柱失稳时对中心土柱的作用力大小，
[0053] τ′Tk为中心土柱计算下滑力时受到的周围土柱的力。
[0054] 可选的，所述利用纤维束模型模拟土柱间的相互作用，并根据连接键断裂数量计算中心土柱与相邻土柱间的应力大小，具体包括：
[0055] 采用纤维束模型模拟土柱间的相互作用，根据纤维束的破坏情况以及土柱间的接触面积，计算土柱破坏过程应力分配的大小；
[0056] 计算公式如下：
[0057]
[0058]
[0059]
[0060] 同理可得，τT1，τT2，τT3，τT5的计算公式，式中n＝0,1,2,3,4,5；
[0061] 式中：Acn表示的是土柱与周围六个的土柱间的接触面积，τT0、τT1、τT2、τT3表示的是i,j土柱对周围土柱的拉应力，τT4、τT5表示的是i,j土柱对周围土柱的压应力，Nn表示i,j土
柱与编号为n的土柱间的完整纤维的数量，K2/K1表示压应力与拉应力间的分配比，hrf表示
中心土柱的侵蚀深度，hrfn表示周围土柱编号为n的土柱的侵蚀深度，H1表示中心土柱坡面
的高程，AH为土柱的截面面积，Lg为相邻土柱间的距离；H1n表示周围土柱编号为n的土柱的
坡面高程，τR为土柱沿最不利滑动面破坏后的残余强度，Hsd为最不利滑动面的深度。
[0062] 一种泥石流坡面物源启动量动态计算系统，所述计算系统包括：
[0063] 坡面划分模块，用于选择待计算的泥石流物源坡面，将所述待计算的泥石流物源坡面划分为规则排列的i行j列的正六边形的土柱；
[0064] 滑动面及不平衡力计算模块，用于根据极限分析上限定理计算每个土柱的最不利滑动面，并计算最不利滑动面上的不平衡力；
[0065] 失稳判断模块，用于根据所述最不利滑动面上的不平衡力分别判断相应土柱是否发生失稳；若判断结果为否，则该土柱停止计算；若判断结果为是，则
[0066] 失稳方向判定模块，用于根据土柱是否失稳，以及中心土柱与其周围六个相邻土柱的重心高程，确定所述中心土柱失稳方向，并获取周围六个相邻土柱对中心土柱的作用
力方式；
[0067] 连接键断裂判断模块，用于利用纤维束模型模拟土柱间的相互作用，并根据连接键断裂数量计算中心土柱与相邻土柱间的应力大小；若判断结果为侧向拉应力的连接键未
全部断裂，则该土柱停止计算；若判断结果为侧向拉应力的连接键全部断裂，则
[0068] 流态化判断模块，用于根据流态化准则判断土柱是否流态化，若判断结果为否，则该土柱停止计算；若判断结果为是，则
[0069] 记录模块，用于记录该土柱启动，记录坡面物源泥石流启动量。
[0070] 可选的，所述滑动面及不平衡力计算模块具体包括：
[0071] 最不利滑动面深度和角度计算单元，用于根据极限分析上限定理，计算最不利滑动面的深度Hsd和角度α+β：
[0072] 总的外功率为：
[0073] 滑动面上的内能耗散为：
[0074] 土柱安全系数K可以表示为：滑动面对应着最小的安全系数：
[0075] 最不利滑动面上的不平衡力计算单元，用于计算最不利滑动面上的不平衡力，具体包括：
[0076] 下滑力计算子单元，用于根据计算最不利滑动面的下滑力；
[0077] 抗滑力计算子单元，用于根据计算最不利滑
动面的抗滑力；
[0078] 最不利滑动面上的不平衡力计算子单元，用于利用所述最不利滑动面的下滑力和所述最不利滑动面的抗滑力计算所述最不利滑动面上的不平衡力；
[0079] 其中： τT为周围相邻6个土柱对中心土柱的应力总和，τ′Tn为所述周围相邻土柱对所述中心土柱的应力；
[0080] Lg表示土柱间的间距，c为待计算的泥石流物源根土复合体的粘聚力，Hsd为最不利滑动面的深度；θ为所述待计算泥石流物源坡面土体含水率，ρw为水的密度， 为所述待计
算泥石流物源坡面根土复合体的孔隙率，ρr为所述待计算泥石流物源坡面根土复合体的密
度，g为重力加速度，α为基岩的角度，β为最不利滑动面同基岩的角度，χ为与水土特征曲线
有关的参数，h为毛细管压力水头，γ为待计算的泥石流物源坡面根土复合体的内摩擦角，
τrf为待计算泥石流物源坡面径流侵蚀力，γw为水的重度，hs为待计算泥石流物源坡面所在
区域的潜水位水头，J为待计算泥石流物源坡面所在区域的潜水位水力坡度，H1为所述中心
土柱的坡面高程，H2为所述中心土柱的基岩高程，Hd为待计算泥石流物源坡面所在区域的地
下水位。
[0081] 可选的，所述待计算的泥石流物源根土复合体的粘聚力c与所述待计算的泥石流物源坡面根土复合体的内摩擦角γ服从威布尔分布：
[0082]
[0083] 所述待计算泥石流物源坡面土体含水率θ、所述待计算泥石流物源坡面根土复合体的孔隙率 与待计算泥石流物源坡面根土复合体的密度ρ服从正态分布：
[0084]2
[0085] 式中：μθ, μρ分别表示含水率、孔隙率和密度服从的正态分布的期望，σθ,2
σρ分别表示含水率、孔隙率和密度服从的正态分布的方差，期望和方差根据坡面取得的实
测值确定。
[0086] 可选的，所述失稳方向判定模块具体包括：
[0087] 中心土柱位置获取单元，用于获取所述中心土柱的位置(i,j)，并确定所述中心土柱的重心坐标为(xi，j,yi,j,zi,j)；
[0088] 周围相邻土柱的重心坐标计算单元，用于根据位置坐标和几何原理得到所述周围相邻土柱的重心坐标(xk,yk,zk)(k＝0,1,2,3,4,5,6)为
[0089]
[0090] 作用向量计算单元，用于根据所述中心土柱的重心坐标和周围相邻土柱的重心坐标得到周围相邻土柱相对中心土柱的向量(x'k,y'k,z'k)，
Lg为相邻土柱间的距离；
[0091] 周围相邻土柱中，两相邻土柱的和向量为(x′k,k+1,y′k，k+1,z′k,k+1)
[0092]
[0093] 所述和向量投影在二维平面内的向量为Dk,k+1(k＝0,1,2,3,4,5)；
[0094] 所述中心土柱运动方向在二维平面投影向量(x″i,j,y″i,j)；
[0095] 土柱运动的方向向量确定单元，用于根据最小作用量原理，取重心最低min(z′k,k+1)(k＝0,1,2,3,4,5)的和向量为土柱运动的方向向量；
[0096] k取min(z′k,k+1)(k＝0,1,2,3,4,5)时的值，所述中心土柱在运动方向在二维平面的投影向量为(x″i,j,y″i,j)：x″i,j＝x′k,k+1、y″i,j＝y′k,k+1；
[0097] 周围土柱运动方向在二维平面投影向量(x″k,y″k)；
[0098] 相邻土柱间的夹角θk为
[0099] 当相邻土柱的运动方向向量的夹角为锐角时，土柱间的作用力为压应力，当相邻土柱间的运动方向向量为钝角的时候，土柱间的作用力为拉应力，
[0100] 若 则：τ′Tk＝‑τTk；
[0101] 若 则：τ′Tk＝0；
[0102] 若 则：τ′Tk＝τTk；
[0103] τTk为周围土柱失稳时对中心土柱的作用力大小，
[0104] τ′Tk为中心土柱计算下滑力时受到的周围土柱的力。
[0105] 可选的，采用纤维束模型模拟土柱间的相互作用，根据纤维束的破坏情况以及土柱间的接触面积，计算土柱破坏过程应力分配的大小；
[0106] 计算公式如下：
[0107]
[0108]
[0109]
[0110] 同理可得，τT1，τT2，τT3，τT5的计算公式，式中n＝0,1,2,3,4,5；
[0111] 式中：Acn表示的是土柱与周围六个的土柱间的接触面积，τT0、τT1、τT2、τT3表示的是i,j土柱对周围土柱的拉应力，τT4、τT5表示的是i,j土柱对周围土柱的压应力，Nn表示i,j土
柱与编号为n的土柱间的完整纤维的数量，K2/K1表示压应力与拉应力间的分配比，hrf表示
中心土柱的侵蚀深度，hrfn表示周围土柱编号为n的土柱的侵蚀深度，H1表示中心土柱坡面
的高程，AH为土柱的截面面积，Lg为相邻土柱间的距离；H1n表示周围土柱编号为n的土柱的
坡面高程，τR为土柱沿最不利滑动面破坏后的残余强度，Hsd为最不利滑动面的深度。
[0112] 根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0113] 本发明将泥石流启动量计算的理论模型数值化，能对潜在的坡面泥石流的启动流量进行准确计算，有利于解决长期存在的泥石流启动量计算难、精确度低的问题；并且运用
纤维束模型模拟土柱和土柱间的破坏及相互作用过程，能够在泥石流启动前根据纤维束的
破坏情况判定土柱的是否即将发生破坏，能对泥石流的启动进行提前预报分析，将有利于
提前规划防灾措施。
附图说明
[0114] 图1为本发明提供的泥石流坡面物源启动量动态计算方法流程示意图；
[0115] 图2为土柱的编号；
[0116] 图3为土柱的位置关系；
[0117] 图4为土柱间应力传递示意图；
[0118] 图5为本发明提供的泥石流坡面物源启动量动态计算系统示意图。
具体实施方式
[0119] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于
本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他
实施例，都属于本发明保护的范围。
[0120] 本发明的目的是提供一种泥石流坡面物源启动量动态计算方法及系统，能对潜在的坡面泥石流的启动流量进行准确计算，将有利于解决长期存在的泥石流启动量计算难、
精确度低的问题。
[0121] 为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0122] 如图1所示，为本发明的一种泥石流坡面物源启动量动态计算方法流程示意图，该计算方法包括：
[0123] 步骤101：选择待计算的泥石流物源坡面，将所述待计算的泥石流物源坡面划分为规则排列的i行j列的正六边形的土柱。
[0124] 步骤102：根据极限分析上限定理计算每个土柱的最不利滑动面，并计算最不利滑动面上的不平衡力。
[0125] 步骤103：根据所述最不利滑动面上的不平衡力分别判断相应土柱是否发生失稳；若判断结果为否，则该土柱停止计算；若判断结果为是，则执行步骤104。
[0126] 步骤104：根据土柱是否失稳，以及中心土柱与其周围六个相邻土柱的重心高程，确定所述中心土柱失稳方向，并获取周围六个相邻土柱对中心土柱的作用力方式。
[0127] 步骤105：利用纤维束模型模拟土柱间的相互作用，并根据连接键断裂数量计算中心土柱与相邻土柱间的应力大小；若判断结果为侧向拉应力的连接键未全部断裂，则该土
柱停止计算；若判断结果为侧向拉应力的连接键全部断裂，则执行步骤106。
[0128] 步骤106：根据流态化准则判断土柱是否流态化，若判断结果为否，则该土柱停止计算；若判断结果为是，则利用步骤107记录坡面物源泥石流启动量。
[0129] 步骤103‑步骤106中的判断结果为否时，也可将坡面物源泥石流启动量启动记录为0。
[0130] 步骤102具体包括：
[0131] 根据极限分析上限定理，最不利滑动面的深度Hsd和角度α+β：
[0132] 总的外功率为：
[0133] 滑动面上的内能耗散为：
[0134] 土柱安全系数K可以表示为：滑动面对应着最小的安全系数：
[0135] 计算最不利滑动面上的不平衡力，具体包括：
[0136] 根据 计算最不利滑动面的下滑力。
[0137] 根据计算最不利滑动面的抗滑力。
[0138] 利用所述最不利滑动面的下滑力和所述最不利滑动面的抗滑力计算所述最不利滑动面上的不平衡力。
[0139] 其中： τT为周围相邻6个土柱对中心土柱的应力总和，τ′Tn为所述周围相邻土柱对所述中心土柱的应力。
[0140] Lg表示土柱间的间距，c为待计算的泥石流物源根土复合体的粘聚力，Hsd为最不利滑动面的深度；θ为所述待计算泥石流物源坡面土体含水率，ρw为水的密度， 为所述待计
算泥石流物源坡面根土复合体的孔隙率，ρr为所述待计算泥石流物源坡面根土复合体的密
度，g为重力加速度，α为基岩的角度，β为最不利滑动面同基岩的角度，χ为与水土特征曲线
有关的参数，h为毛细管压力水头，γ为待计算的泥石流物源坡面根土复合体的内摩擦角，
τrf为待计算泥石流物源坡面径流侵蚀力，γw为水的重度，hs为待计算泥石流物源坡面所在
区域的潜水位水头，J为待计算泥石流物源坡面所在区域的潜水位水力坡度，H1为所述中心
土柱的坡面高程，H2为所述中心土柱的基岩高程，Hd为待计算泥石流物源坡面所在区域的地
下水位。
[0141] 上述参数中，待计算的泥石流物源根土复合体的粘聚力c与待计算的泥石流物源坡面根土复合体的内摩擦角γ服从威布尔分布：
[0142]
[0143] 待计算泥石流物源坡面土体含水率θ、待计算泥石流物源坡面根土复合体的孔隙率 与待 计算 泥石 流 物源 坡面 根土 复 合体的 密 度ρ服 从正 态分 布 ：
2
[0144] 式中：μθ, μρ分别表示含水率、孔隙率和密度服从的正态分布的期望，σθ,2
σρ分别表示含水率、孔隙率和密度服从的正态分布的方差，期望和方差根据坡面取得的实
测值确定。
[0145] 待计算泥石流物源坡面根土复合体的密度ρ、孔隙率 粘聚力c和内摩擦角γ等，是通过选取待计算泥石流物源坡面根土复合体，通过室内土工试验测量得到的。
[0146] 所述待计算泥石流物源坡面土体含水率θ、水土特征曲线参数λ、待计算泥石流物源坡面所在区域的潜水位水头hs、和潜水位水力坡度J是采用相关理论或经验公式如
Richard方程、Darcy定律等计算得到的。
[0147] 步骤104具体包括：
[0148] 土柱失稳破坏运动的倾向，决定了土柱间的作用力分配方向。根据最小作用量原理，可以根据破坏土柱的重心判定土柱的运动倾向。图2和图3为土柱的编号及位置关系，i,
j分别表示土柱的位置处于i行j列。
[0149] 对周围土柱进行编号，如图2所示。D0,1、D1,2、D2,3、D3,4、D4,5、D5,6表示六个土柱发生失稳破坏的6个倾向。以D5,6为例，若中心土柱i,j沿D5,6方向发生失稳破坏，则沿D5,6方向i,j
土柱前方的周围土柱i,j+2和i+1,j+1土柱受到，i,j土柱的压应力，其他土柱则受到i,j土
柱的拉应力。
[0150] 获取中心土柱的位置(i,j)，并确定中心土柱的重心坐标为(xi，j,yi,j,zi,j)。
[0151] 根据位置坐标和几何原理得到所述周围相邻土柱的重心坐标(xk,yk,zk)(k＝0,1,2,3,4,5,6)为
[0152]
[0153] 根据所述中心土柱的重心坐标和周围相邻土柱的重心坐标得到周围相邻土柱相对中心土柱的向量(x'k,y'k,z'k)，为 Lg为相邻土
柱间的距离。
[0154] 周围相邻土柱中，俩相邻土柱的和向量为(x′k,k+1,y′k,k+1,z′k,k+1)，(k＝0,1,2,3,4,5)
[0155]
[0156] 所述和向量投影在二维平面内的向量为Dk,k+1(k＝0,1,2,3,4,5)。
[0157] 所述中心土柱运动方向在二维平面投影向量(x″i,j,y″i,j)。
[0158] 根据最小作用量原理，取重心最低min(z′k,k+1)(k＝0,1,2,3,4,5)的和向量为土柱运动的方向向量。
[0159] k取min(z′k,k+1)(k＝0,1,2,3,4,5)时的值，所述中心土柱在运动方向在二维平面的投影向量为(x″i,j,y″i,j)：x″i,j＝x′k,k+1、y″i,j＝y′k,k+1。
[0160] 周围土柱运动方向在二维平面投影向量(x″k,y″k)。
[0161] 相邻土柱间的夹角θk为
[0162] 当相邻土柱的运动方向向量的夹角为锐角时，土柱间的作用力为压应力，当相邻土柱间的运动方向向量为钝角的时候，土柱间的作用力为拉应力。
[0163] 若 则：τ′Tk＝‑τTk。
[0164] 若 则：τ′Tk＝0。
[0165] 若 则：τ′Tk＝τTk。
[0166] τTk为周围土柱失稳时对中心土柱的作用力大小。
[0167] τ′Tk为中心土柱计算下滑力时受到的周围土柱的力。
[0168] 在纤维束模型中，最不利滑动面或土柱接触面间的荷载为： 纤维束的强度为： 完整纤维束的数量为：
[0169] 式中：NI为完好纤维束的数量，NF为接触面上纤维束的总数量，fR为纤维束破坏后的残余强度系数(fR＝纤维束断裂后的残余抗拉强度/纤维束断裂前的抗拉强度)，σF为施加
在每根完好纤维束上的力，σm狐x(θ)为一根纤维束所能承受的最大抗拉强度和含水率相关
(通过实验测量)，Ac为接触面的面积，τh(θ)为根‑土复合体的抗拉强度，与含水率有关(通过
实验测量)。
[0170] NI表示的是滑动面或土柱接触面上的完整纤维束的数量，NF表示的是滑动面或土柱接触面上总的纤维束的数量，从上述的公式可以发现，NI的大小受到接触力的影响，NI的
大小直接反映了土柱基底的破坏情况。当NI＝NF时，表示土柱破坏面上纤维束未发生断裂，
即土柱完整，还未发生任何损伤破坏；当NINI和NF的比值达到破坏阈值时(根据具体的坡面设置阈值)，即可认为土柱已经发生破坏，即
将物源启动。通过记录NI的值，我们可以预知土柱的破坏程度，即可以预测坡面潜在的易发
生破坏的区域，在坡面发生破坏前采取相应的工程措施，从而达到防灾减灾的目的。
[0171] 利用纤维束模型模拟土柱间的相互作用，根据纤维束模型判断相邻土柱接触面上总的纤维数量与完整纤维束数量的比值是否大于设定阈值；若判断结果为否，则侧向拉应
力的连接键未全部断裂；若判断结果为是，则侧向拉应力的连接键全部断裂。并根据纤维束
的破坏数量情况以及土柱间的接触面积，计算土柱破坏过程应力分配的大小。
[0172] 当坡面的土柱单元沿最不利滑动面发生失稳后，土柱在根系、胶结、孔隙水等的作用下，通过连接键向沿着土柱破坏运动方向的前方土柱分配压应力，向沿着土柱破坏运动
方向的后方和侧方土柱分配拉应力。
[0173] 应力传递过程如图4所示，计算公式如下：
[0174]
[0175]
[0176]
[0177] 同理可得，τT1，τT2，τT3，τT5的计算公式，式中n＝0,1,2,3,4,5。
[0178] 式中：Acn表示的是土柱与周围六个的土柱间的接触面积，τT0、τT1、τT2、τT3表示的是i,j土柱对周围土柱的拉应力，τT4、τT5表示的是i,j土柱对周围土柱的压应力，Nn表示i,j土
柱与编号为n的土柱间的完整纤维的数量，K2/K1表示压应力与拉应力间的分配比，hrf表示
中心土柱的侵蚀深度(采用Iverson等提出侵蚀模型进行计算)，hrfn表示周围土柱编号为n
的土柱的侵蚀深度，H1表示中心土柱坡面的高程，AH为土柱的截面面积，Lg为相邻土柱间的
距离。H1n表示周围土柱编号为n的土柱的坡面高程，τR为土柱沿最不利滑动面破坏后的残余
强度。
[0179] 土柱流态化的判断依据为： 其中τc为中心土柱受到的抗压强度，τh(θ)为中心土柱同含水率有关的抗剪强度；γ为待计算的泥石流物源坡面根
土复合体的内摩擦角。
[0180] 如图5所示，为本发明的一种泥石流坡面物源启动量动态计算系统示意图，该计算系统包括：
[0181] 坡面划分模块201用于选择待计算的泥石流物源坡面，将所述待计算的泥石流物源坡面划分为规则排列的i行j列的正六边形的土柱。
[0182] 滑动面及不平衡力计算模块202用于根据极限分析上限定理计算每个土柱的最不利滑动面，并计算最不利滑动面上的不平衡力。
[0183] 失稳判断模块203用于根据所述最不利滑动面上的不平衡力分别判断相应土柱是否发生失稳；若判断结果为否，则该土柱停止计算，连接停止模块208；若判断结果为是，则
失稳判断模块203连接失稳方向判定模块204。
[0184] 失稳方向判定模块204用于根据土柱是否失稳，以及中心土柱与其周围六个相邻土柱的重心高程，确定所述中心土柱失稳方向，并获取周围六个相邻土柱对中心土柱的作
用力方式。
[0185] 连接键断裂判断模块205用于利用纤维束模型模拟土柱间的相互作用，并根据连接键断裂数量计算中心土柱与相邻土柱间的应力大小；若判断结果为侧向拉应力的连接键
未全部断裂，则该土柱停止计算，连接停止模块208；若判断结果为侧向拉应力的连接键全
部断裂，则连接键断裂判断模块205连接流态化判断模块206。
[0186] 流态化判断模块206用于根据流态化准则判断土柱是否流态化，若判断结果为否，则该土柱停止计算，连接停止模块208；若判断结果为是，则利用记录模块207记录该土柱启
动，记录坡面物源泥石流启动量。
[0187] 上述模块203‑206可还分别与记录模块207连接，在判断结果为否时，将坡面物源泥石流启动量启动记录为0。
[0188] 对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0189] 本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据
本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不
应理解为对本发明的限制。
