黄土沟壑区沟坡稳定分析及失稳概率计算

黄土沟壑区沟坡稳定分析及失稳概率计算

重力侵蚀是指在重力的作用下，边坡上的岩石和土壤失去平衡，移动时不可避免地会发生位移变化的过程。重力侵蚀包括腹泻、滑动、坍塌、坍塌和泥石流等。在我国的黄土高原地区，重力侵蚀的类型多种多样。根据刘秉正对渭北地区的调查，导致腹泻滑动的频率最高，约为56.5%，但侵蚀的最小频率约为3.0%。虽然滑动和崩塌的发生频率不高，但侵蚀超过96%。因此，滑坡和崩塌是主要的重力侵蚀类型。本文以物理机为基础，建立了这种重力侵蚀的计算模型。重力侵蚀不同于普通工程意义上的滑坡失稳,主要有如下两方面的原因.首先,黄土丘陵沟壑区的沟坡一般是由于坡面下部切沟在水流作用下不断下切,切沟之间形成了孤立沟坡,而沟坡在降雨等条件下,由于重力作用产生崩塌.另一方面,沟坡下部在雨季形成行洪沟道,沟坡本身将受到水流的强烈淘刷而加大直立面高度,从而对沟坡的稳定性产生严重影响.因此,黄土沟壑区的重力侵蚀是多因素作用下的坡面失稳问题,水流的切割和降雨影响必须作为主要因素加以考虑.由于重力侵蚀的复杂性,以及缺乏长期系统的观测资料,对重力侵蚀的研究在土壤侵蚀领域中一直是一个较薄弱的环节.近年来,也有学者运用各种方法对重力侵蚀进行研究.李容全等利用遥感影像提取沟谷形态,对米脂地区重力侵蚀的平均侵蚀模数和侵蚀强度进行了计算.朱同新,陈永宗等利用模糊聚类分析的方法将晋西北的重力侵蚀分为5个区,并对各区的重力侵蚀特征进行了分析.王军等通过扩展地理侵蚀面的概念,利用GIS技术,提出了黄土高原地区重力侵蚀相对强度的估算式,并对相对强度的空间分布进行了探讨.上述研究者大多是从地理地貌学或者利用GIS的空间分析功能方面对重力侵蚀特征进行研究,而重力侵蚀的发生最根本原因在于其力学机理.本文就是在这方面进行探讨,提出基于物理机理的重力侵蚀力学模型.1重力侵蚀模型的建立1.1产流模型和数值模型重力侵蚀通常发生在黄土坡面的沟谷区,其形态十分复杂.本文是针对数字流域模型中关于流域单元划分的第四级,即坡面上建立重力侵蚀模型.因此,根据数字流域模型平台提取的坡面信息,将重力侵蚀发生的单位宽度坡面概化如图1所示.图中将坡面土体分为表层土和中层土,其入渗过程由相应的产流模型计算.沟坡受力状况如下:(1)土体重力Wt.本模型中土体重力主要考虑土体自重以及降雨过程中由于雨水入渗产生的重力增量;(2)裂缝水压力T.由于黄土垂直节理较为发育,坡面常形成垂直裂缝.降雨条件下,雨水灌入将形成一定的水平水压力,对土体稳定造成不利影响;(3)裂隙面上的抗滑力FR.抗滑力与土体的粘聚力和摩擦角有关.降雨条件下,随含水量变化;(4)坡脚处的水流切割力τ.该力将直接导致坡角的横向后退,继而使直立面高度增加,对坡面稳定产生不利影响.沟坡重力侵蚀是在上述各力的作用下,达到临界平衡条件导致的失稳破坏.按照失稳破坏的实际发生过程,需同时模拟水流切割以及抗剪力减小导致的稳定性降低,然后运用模糊概率分析的方法,把安全系数转化为失稳概率进行预报,具体流程如图2所示.1.2模型的构建1.2.1b.地形地貌条件黄土坡脚下形成的沟道是雨季主要的行洪通道.由于其坡降一般较大,洪峰较为集中,因此水流对坡脚将形成强烈的侧向淘刷.黄土沟坡颗粒粒径在0.05—0.005mm之间的粉粒占65%左右,而粒径小于0.001mm的粘性颗粒约占6%,因此,物理特性可按具有一定粘结力的粘土考虑.唐存本将重力、拖曳力、上举力、粘结力统一考虑,得出新淤粘土的起动切应力公式,本模型的起动切应力应用此公式:τc=66.8×102×d+3.67×10−6d(1)τc=66.8×102×d+3.67×10-6d(1)式中τc为起动切应力(N/m2);d为粒径(m).在给定的Δt时间内,洪水持续对沟坡进行侧向冲刷,冲刷导致的横向后退速度与水流切应力τ及上述的土体起动切应力τc有关,同时还与土体本身的理化性质有关.本文采用Osman等根据室内模型试验得到的土体单位时间侧向冲刷距离:ΔB=Cl×(τ−τc)×e−1.3τcγs(2)ΔB=Cl×(τ-τc)×e-1.3τcγs(2)式中ΔΒ为土体单位时间受水流侧向冲刷而后退的距离(m);τ为水流切应力(N/m2);τc为土体起动切应力(N/m2);Cl为与土体理化特性有关的系数,根据Osman的试验资料,可取Cl=3.64×10-4.由图1所示的几何关系可知,坡脚由于水流侧向冲刷而后退ΔΒ后,沟坡将相应产生直立高度,其转折点之上的沟坡高度为:H1=H−ΔBtaniΗ1=Η-ΔΒtani式中H1为直立面转折点上的沟坡高度(m);H为沟坡高度(m);i为沟坡自然坡角度;ΔΒ同前.当沟坡发生垮塌时,破坏面与水平面的夹角为:β=0.5×{tan−1[(HH1)(1−k2)tan(i)]+ϕ}(3)β=0.5×{tan-1[(ΗΗ1)(1-k2)tan(i)]+ϕ}(3)式中k为黄土坡面中较大垂直节理或裂隙深度与沟坡高度H的比值,可根据地质调查确定,无资料时可取0.3;ϕ为摩擦角;β为垮塌面与水平面的夹角;其余各量同上.1.2.2抗滑性能的确定(1)沟坡滑动强度确定降雨期间,沟坡的下滑力主要由坡面垂直节理中水压力、沟坡土体重力、入渗雨水重力等在破坏面上的分力组成.1直节理裂隙深度的计算由于黄土普遍具有较为发育的垂直节理,在各中营力作用下,黄土坡面,特别是坡度较陡的沟坡部分常沿某一垂直节理产生具有一定深度的裂隙.在降雨期间裂隙中充满雨水,对深度较大的裂隙中的水压力,是不可忽略的,将构成沟坡下滑力的一部分.该水压力在沿破坏面的分力可由下式计算:T=12γH2t(4)Τ=12γΗt2(4)式中γ为入渗雨水容重(kN/m3);Ht为裂隙深度(m);T为入渗雨水压力(kN/m).2土体重力wt的计算由图1概化沟坡的几何关系可知,沟坡土体的重力可表示为:Wt=γwm2(H2−H2ttanβ−H21tani)(5)Wt=γwm2(Η2-Ηt2tanβ-Η12tani)(5)式中Wt为可能失稳的土体重力(kN/m);γwm为相应于某一土体含水量w时的土体容重(kN/m3);其余各量同前.3沟坡智能分力组成下滑力由上述裂隙中水压力和土体重力沿失稳破坏面的分力组成,可由下式表达:FD=Wtsinβ+Tcosβ(6)FD=Wtsinβ+Τcosβ(6)式中FD为沟坡下滑力(kN/m);其余各量同前.(2)非饱和土抗剪强度的确定抗滑力的确定较为复杂.由于黄土高原地区的降雨一般都集中在汛期的几个月时间里有限的几场暴雨中,同时由于地下水位较低,因此,大部分时间沟坡土体均处于非饱和状态.非饱和土体的本构关系与强度特征与饱和土有较大不同.主要是非饱和土体不仅要满足土体本身的应力应变关系,同时还受到土体含水量的较大影响.在含水量较低时,由于负空隙水压力的存在,形成基质吸力,从效果上看,相当于增加了附加粘聚力,增强了土体的抗剪强度.但随着降雨入渗,土体的含水量增大,则附加粘聚力急剧降低,从而导致抗剪强度的减小,当抗剪强度不足以抵抗下滑力时,就可能发生沟坡的滑动破坏.对于基质吸力随含水量变化的研究,是非饱和土力学中的重要课题.但由于基质吸力与土壤种类、土壤特性、应力条件、初始条件等具有复杂的关系,因此至今还未形成统一的结果.对黄土基质吸力和强度的研究比较活跃.胡在强、沈珠江等通过轴平移法量测吸力技术,研究了人工制备的结构性黄土吸力随含水量及饱和度的变化关系.张伯平等对杨陵渭北黄土的结构强度与含水量的关系进行了试验研究.党进谦等分析了非饱和黄土强度的组成和来源,证明了黄土结构强度与初始含水量具有幂函数关系,并给出了非饱和黄土抗剪强度的拟合公式.目前得到广泛承认的非饱和土抗剪强度的理论表达式,由Fredlund等于1978年提出,如下式所示:τ=c′+(σ−ua)tanϕ′+(ua−uw)tanϕb(7)τ=c′+(σ-ua)tanϕ′+(ua-uw)tanϕb(7)式中τ为非饱和土的抗剪强度;c′和ϕ′为饱和土体的有效粘聚力与摩擦角;σ为总应力;ua为孔隙气压力;uw为孔隙水压力;ua-uw即为基质吸力;ϕb为抗剪强度随基质吸力增加的斜率,其本身又是基质吸力的函数.由上式可知,要确定非饱和土的抗剪强度,必须以确定基质吸力及剪切力增量的斜率随含水量的变化为前提,而他们本身又互为函数,因此按其定义式,目前较难应用于实际计算.另一种简化的处理方法是,将非饱和土抗剪强度中的粘聚力分为饱和粘聚力与附加粘聚力,并通过试验得到不同土体的附加粘聚力随含水量的变化关系.卢肇钧指出土体的内摩擦角随含水量的变化一般不超过±2°,因此,在计算中对于同一土体,其内摩擦角可以按常数考虑.由此,非饱和土抗剪强度可近似写为:τ=c+σtanϕ=c′+τ′+σtanϕ(8)τ=c+σtanϕ=c′+τ′+σtanϕ(8)式中τ为非饱和土的抗剪强度(kPa);c′为相应于饱和土体的粘聚力(kPa);τ′为附加粘聚力;其余各量同前.党进谦等对渭北台塬地的Q2和Q3非饱和黄土,以及关中地区的十个取土点的黄土的结构强度、抗剪强度、基质吸力等与含水量的关系进行了较为深入的分析.实验数据表明,附加粘聚力与含水量之间存在良好的幂函数关系:τ′=awb(9)τ′=awb(9)式中τ′为非饱和黄土的附加粘聚力(kPa);w为土体含水量(%);a、b为幂函数的系数和指数,通过试验得到.根据党进谦的试验资料,对于Q2,Q3黄土,经回归得到如下拟合图像与关系:对于不同地区的黄土,其附加粘聚力随含水量变化的系数及指数有一定变化,应根据相关试验拟合确定.确定黄土抗剪强度参数随含水量变化的关系后,即可由图1所示的几何关系,确定沟坡滑动面上所受的抗滑力,如下式所示:FR=cL+Ntanϕ=(H−Ht)csinβ+Wtcosβtanϕ(10)FR=cL+Νtanϕ=(Η-Ηt)csinβ+Wtcosβtanϕ(10)式中FR为滑动面上的抗滑力(kN/m);N为作用在滑动面上的法向力(kN/m);c由式(8)、式(9)确定;其余各量同前.2.2.3抗滑力与下滑力按照传统的极限平衡法,定义安全系数为抗滑力与下滑力的比值,即:Fs=FRFD(11)Fs=FRFD(11)式中FR,FD分别为按式(10)和式(6)计算出的沟坡抗滑力与下滑力;Fs为沟坡抗滑安全系数.1.2.4预报条件分析使用上述方法确定沟坡的安全系数Fs后,仍然不能进行正确的沟坡重力侵蚀模拟,主要是由于以下原因:(1)安全系数只是安全性的一种特征,用于工程安全评价时是可行的,但用于沟坡失稳分析时,不能直接作为沟坡失稳的预报条件;(2)从严格的意义上看,安全系数的各制约因素,例如粘聚力、摩擦角等即使对于同一土体,也不是完全确定的值,而是随机变量;(3)安全系数概念本身就具有一定的模糊性,例如安全系数大于1并不等价于沟坡不失稳这一命题.不确定性包括两个方面,一方面是随机性,另一方面是模糊性.在本问题中,如果定义“该安全系数表示沟坡失稳”为事件S,则该事件S为模糊事件,因为对于某一确切的Fs,并不能确定沟坡是否确实失稳.同时,陈立宏,陈祖煜等通过收集多个工程的长系列抗剪强度试验资料发现,各抗剪强度指标均可接受正态分布和对数正态分布.因此,沟坡失稳问题可以归结为模糊事件的概率预报问题.(1)偏小型模糊分布安全系数模糊程度可用模糊数学中的隶属度函数进行刻画.设安全系数Fs在某个区间(如(0,1.5])内连续变化,则该区间构成问题的一个论域F.于是定义:A:F→,Fs|→A(Fs)(12)A:F→,Fs|→A(Fs)(12)则:A为从论域F到[0,1]上的映射,构成模糊集,且其意义表示Fs在论域中变化时,沟坡失稳的可能性;A(Fs)为定义在论域F上的隶属度;当Fs为连续型变量时,A(Fs)为隶属度函数.隶属度函数具有多种形式,在本问题中,由于失稳隶属度随着安全系数的增大而趋于0,可以选择的类型为偏小型模糊分布.其子类型包括:降半正态分布、将半Cauchy(柯西)分布、降岭型分布、Z型分布等.假设在(0,1.5]的区间内隶属度函数满足:当0<Fs<0.8有A(Fs)=1;当0.8<Fs<1.5有0<A(Fs)<1;且A(Fs)为单调减函数;当Fs>1.5有A(Fs)=0;则通过对比计算,得到上述四种分布在该区间的形态,如图4所示.需要说明的是,对于Z分布,经过数值试验,需要将原分布中的2次方替换为4次方,即如下式所示:A(Fs)=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪11−12(Fs−F1F2−F1)4−12(Fs−F3F3−F2)400<Fs<F1F1<Fs<F2F2<Fs<F3Fs>F3(13)A(Fs)={10<Fs<F11-12(Fs-F1F2-F1)4F1<Fs<F2-12(Fs-F3F3-F2)4F2<Fs<F30Fs>F3(13)式中F1、F2、F3分别为Fs在论域上相对于隶属函数值为1,中值和0值的分界点.由图中可以看出,方次修改为4次的Z分布,当安全系数Fs在论域中减小时,其失稳隶属函数值趋近于1的速度最快,而当安全系数Fs在论域中增大时,其失稳隶属函数趋近于0的也较大,这种函数形态比较符合实际中当安全系数减小到一定程度,失稳可能性急剧增加,当安全系数增大到一定程度,失稳可能性急剧降低的经验,因此,在上述安全系数变化区间内,方次为4的Z分布更适于作为失稳隶属函数.(2)模糊集合空间的概率测度沟坡失稳的随机性表现在,其各种影响因素都具有随机性,这些随机因素又相互消涨,使沟坡稳定性总体上也表现出随机性.罗文强等将一次二阶矩方法应用于斜坡稳定性概率评价中,并讨论了安全系数与破坏概率的关系;徐卫亚等利用蒙特卡罗方法对破坏概率进行了分析,并对水布垭坝址基岩的稳定性进行了计算;罗文强等对矩估计、蒙特卡罗法、函数连分式等几种斜坡可靠度模型进行了对比分析,指出这些方法在精度和难易程度上各有特点.实际上,如果将普通集合空间扩展到模糊集合空间,则可定义模糊事件的概率测度如下式所示:Pz(A)=∫ΩA(Fs)dP=∫ΩA(Fs)f(Fs)dFs(14)Ρz(A)=∫ΩA(Fs)dΡ=∫ΩA(Fs)f(Fs)dFs(14)式中Pz(A)为模糊事件A的概率;Ω为概率空间;A(Fs)为安全系数Fs的隶属度函数,由式(13)确定;f(Fs)为安全系数