基于降雨入渗的Newmark模型改进及地震滑坡危险性预测研究

雷 真,李林锐,隆交凤,陈敬男,杨 洋

(云南大学 建筑与规划学院,云南 昆明 650500)

0 引言

滑坡是地震与降雨诱发的一种主要次生灾害,造成的经济损失往往比地震或强降雨造成的直接损失大得多。导致滑坡发生的因素有多种,其中降雨与地震两种因素尤为重要。降雨渗透将导致斜坡土体含水量与重度增加,下滑力增大,抗剪强度降低。尤其在土体与岩体的交界面上,由于土、岩交界地层形成的年代普遍较晚。上部覆盖的土层孔隙较大,也较松散,下部岩层又表现出较强的力学性质,在降水渗透的情况下,更容易导致滑坡失稳。Green-Ampt(GA)降雨入渗模型是目前降雨滑坡分析中应用最为广泛的模型之一[1]。传统模型计算简单、易于求解,但是该模型求解的是水平面,而未考虑斜面。目前基于该模型,大量学者对其进行改进修正。针对初始含水率因素,王千等[2]提出了初始含水率为线性分布时的GA入渗模型,Chen等[3]提出了初始含水率均匀的斜面降雨入渗GA模型。彭振阳等[4]基于分层假设的GA模型,引入湿润层等效导水率改进了传统GA入渗模型。刘卫涛等[5]考虑土体非饱和特性,对传统GA模型进行修正,并应用到斜坡稳定性研究中。李秀珍等[6]结合Mein-Larson入渗模型进行了饱和土、非饱和土与近似饱和土三种情况下的斜坡稳定性分析,对Green-Ampt入渗模型应用与饱和土与近似非饱和土斜坡稳定性模型进行了改进。

Newmark累积位移模型[7]是Newmark于1965年首次提出且随着GIS技术的不断创新和发展,该方法被国内外学者大量的应用于地震滑坡危险性评估研究中[8-12]。Newmark累积位移回归方程可基于强震记录数据进行经验拟合,Jibson[8]根据世界范围内30次地震的2 270条强震记录得到的数据建立了多种位移回归方程;徐光兴等[9]采用Newmark方法基于强震动记录提出斜坡永久位移预测模型,验证了采取CAR临界加速度比为参数的模型更加适用。许多学者采用Newmark方法对“5·12”大地震受灾区域进行了地震滑坡危险性评估分析与区划图的编制[10-11]。刘甲美等[12]基于Newmark模型开展了九寨沟地震滑坡快速评估,并得出该次地震滑坡主要以浅表型碎屑流和小规模崩塌居多的结论。金凯平等[13]提出对Newmark累积位移的修正,并将该模型应用于芦山地震进行地震滑坡危险性区划。

然而目前采用Newmark模型进行地震滑坡危险预测中考虑降雨因素的研究较少,且降雨和地震耦合作用必将提高滑坡发生的概率,因此分析两因素同时诱导滑坡发生的情况是非常有必要的。本文基于GA降雨入渗模型提出了修正的Newmark累积位移模型,对斜坡安全系数Fs进行推导。以云南省鲁甸县某区域为例开展降雨和地震耦合作用下的滑坡危险性预测,并对预测结果进行分析,为今后的滑坡危险性预测提供一定的理论依据。

1 GA模型修正及地下水位确定

1.1 GA模型修正

降雨入渗后,地下水位可以通过一个水文模型来描述,如图1所示。假设该水文模型的初始含水率均匀分布,图中α为斜坡坡脚,q为降雨强度,则可通过Green-Ampt降雨入渗模型并简化为一维入渗问题推导湿润锋深度h。

降雨累积入渗量与湿润锋深度可通过质量守恒理论表示为:

(1)

式中:F为降雨入渗总量;θs为土体饱和含水率;θi为土体初始含水率;h为湿润锋深度。

在降雨强度q小于土体饱和导水系数KS时,可认为降雨全部被土体吸收,湿润锋深度h可表示为:

(2)

如果降雨强度q大于土体饱和渗透系数KS,随着降雨时间的增长,地面将在某个时刻开始积水,假设该临界时刻为降雨历时t=tp且湿润锋深度以上土体均为饱和土体,则在积水前入渗深度可由式 (2)计算;令Sf为湿润峰深度h处的基质吸力,则积水后垂直于斜面的降雨入渗深度可表示为:

(3)

Sf可由王红闪等[14]提出的方法计算,即:

(4)

式中:S为土壤吸湿率。因此降雨入渗深度h随时间t增加可表示为:

(5)

根据式(5)可以得到在临界时刻tp时,由于降雨入渗深度为关于时间t的连续型函数,将tp带入式(5)得到。

(6)

将式(3)代入整理可以得到临界时刻tp的关系式。

(7)

1.2 地下水位的确定

随着降雨入渗过程的发展,饱和土体内地下水将沿着渗流管网平行于坡面进行流动,因此湿润锋深度也随之发生变化,此过程可假设平行坡面的水流由一个体积为V的容器控制,则对于该体积,假定土体性质是恒定的,则在时间t流出体积V的水流量可以用高度变化Δh(图1)来表示:

(8)

式中:x为元素单位大小;s是元素单位间的距离。

根据Beven提出的饱和带内任何点的水力梯度等于模型坡度的假设,则V可表示为:

V=KShstsinα

(9)

则修正后的入渗深度:

Hm=h±Δh

(10)

本文基于ArcGIS平台采用D8单流向算法按上述式子计算Hm。D8流向算法即8个方向建模,是ArcGIS平台水文分析重要的工具,其工作原理如图2所示。本文假设雨水入渗至土体与岩石交界面处不再下渗,最后将修正后的入渗深度视为地下水位深度Hw,即Hw=Hm。

图2 D8算法工作原理Fig.2 Working principle of D8 algorithm

2 Newmark模型的修正

Newmark累积边坡位移模型是Newmark于1965年提出的一种极限平衡理论模型,虽然该模型未考虑竖向地震作用,但仍在地震滑坡危险性预测分析中被广大学者使用。其原理是将边坡视为沿斜面滑动的刚块,当外力达到使边坡由稳定到失稳的最小力时,边坡开始滑动,在地震发生时可将地震引起的外力用mac表示且方向平行于坡面。在降雨期间,随着雨水入渗,土体含水量逐渐增加,斜坡抗剪强度降低,孔隙水压力不断增大,致使滑坡发生的概率增加。本文假设在地震发生时孔隙水压力不随着地震作用发生变化,不考虑地震所产生的超孔隙水压力,并将Green-Ampt模型与Newmark模型相结合,对Newmark模型进行改进(图3)。

图3 改进的Newmark模型Fig.3 Improved Newmark model

在无降雨期间,坡体可视为干燥条件不考虑孔隙水压力的作用,则静态安全系数FS可表示为:

(11)

(12)

在降雨期间,随着降雨时间的增加,孔隙水压力不断增大,此时安全系数FS采用如下式计算:

FZ=csatb+(P-U)tanφsatcosα

(13)

FX=Psinα

(14)

P=b[γ(H-Hw)+Hwγsat]

(15)

U=γwbHw

(16)

(17)

地震作用下,当坡体达到极限平衡时,即FS=1,则临界加速度ac可表示为:

ac=(FS-1)gsinα

(18)

式中:P为坡体重度(kN);H为坡体厚度(m);U为孔隙水压力(kN);c和csat分别为土体的天然和饱和黏聚力(MPa);φ和φsat分别为土体的天然和饱和内摩擦角(°);m为坡体质量;g为重力加速度;ac为致使坡体下滑失稳的临界加速度;γ、γsat和γw分别表示土体天然重度、土体饱和重度与水的重度(kN/m3);b为坡体长度(m)。

目前计算滑坡累积位移DN的模型有很多,代表性的模型主要分为三种,第一种是积分模型:

DN=∬[a(t)-ac]dtdt

(19)

式中:a(t)为加速度时程,由于中强震动记录比较缺乏且大量精确积分难度较高,因此该模型应用较为不便;第二种模型则是以临界加速度比(CAR)为参数计算的模型;第三种则是以Arias强度与CAR为参数的计算模型。通过事实分析,Arias强度虽能较好地反映出地震动特性,但其获得的过程复杂且精度未见提高,因此不建议使用[9]。

lg(DN)=0.546MW-0.425CAR-1.086±0.263

(20)

3 模型应用

3.1 试验区基本地质特征及计算参数

本文选取鲁甸县某个区域作为试验区进行该方法的滑坡预测应用。鲁甸县地处云贵高原西北部,滇东北高原南部(103°09′～103°40′E,26°59′～27°32′N)[15],为昭通市一县城。南北向构造与东北向构造共同组成了鲁甸县的主干构造。其中,南北向构造的主要断层与褶皱有野牛塘断层、五里牌断层、韦家渡背斜、新街子向斜等,东北向构造主要为骡马口构造带[16]。鲁甸县区岩性主要为灰岩、白云岩、砂岩、泥岩、页岩和玄武岩,其中玄武岩大片分布。泥岩与页岩岩性较为软弱,加之玄武岩风化严重,导致研究区地质环境整体脆弱,为滑坡发育提供基础[17]。根据当地地质特征及降雨情况资料,该试验区相关计算参数选取列于表1。

通过该研究区Dem(30 m)数据生成坡度图,图4所示。根据鲁甸县以往的地震记录,本次模拟取MW=6.1。根据鲁甸地震时,美国地质勘探局(USGS)所公布的PGA(PGA，地震动峰值加速度)数据如图5所示。由于地震诱发的滑坡大多数为浅层滑坡,根据前人经验取滑坡厚度H=3 m[10-12]。

图4 坡度图[单位:(°)]Fig.4 Slope map [Unit:(°)]

图5 PGA图(单位:g)Fig.5 PGA map (Unit:g)

3.2 滑坡位移计算及危险区划

通过假设该试验区在均衡降雨强度q作用下,由表1参数可知假设的降雨强度q大于土体饱和渗透系数KS,则存在临界降雨时刻tp使得地面开始积水,通过1.1节公式可推导tp:

表1 研究区参数表Table 1 Parameters of study area

(21)

则通过ArcGIS平台计算得出研究区开始积水的临界时刻tp(图6),研究区中白色区域为坡度小于10°的区域,依据经验[13],该部分区域不发生滑坡,因此不参与计算。

图6 临界时刻tp图(单位:h)Fig.6 Critical time tp map (Unit:h)

为方便计算,本文取该研究区tp数据的中值tp=5 h,分别计算降雨历时5 h(此时地面还未积水)与降雨历时10 h地下水位深度Hw(图7)。

图7 地下水位深度图(单位:m)Fig.7 Depth of groundwater level (Unit:m)

通过计算分别得到无降雨时期(t=0 h)、降雨历时t=5 h(地面无积水)以及降雨历时t=10 h(地面积水)的地震滑坡累计位移Dn。Dn值的大小与滑坡发生的危险性呈正相关的趋势,参考前人经验并根据以往实际滑坡的观测数据可以发现,在安全区域或少量发生滑坡的区域,其累积位移计算值往往小于0.5 cm,危险程度较低;而发生大规模地震滑坡灾害的区域,其累积位移值均处于2 cm以上[12,18],所以将累积位移值分为三类“DN<0.5 cm(危险程度低)、0.5 cm2 cm(危险程度高)”。最终结果如图8所示。

图8 滑坡预测区划图Fig.8 Landslide prediction zoning map

3.3 滑坡危险区分布特征及结果分析

依据滑坡预测区划图可知,滑坡高危险区域主要集中在PGA与坡度值高的地区,表明地震滑坡危险程度与坡度和PGA具有较高的相关性。随着降雨时间的增加,危险程度低的部分区域逐步转变为危险程度中的区域,部分中危险程度区域转变为高危险区域。在降雨期间,随着降雨时间的增加(t=0 h→10 h),区域危险性程度不断上升,危险面积占比同步提高。与未降雨情况相比,降雨历时5 h、10 h情况下地震滑坡低危险区域占比面积从51%分别下降至35%、33%,下降幅度分别为31%、35%。高危险区域面积从1%分别提高至9%、12%,提高幅度分别为800%、1 100%。中危险区域面积占比从48%分别增加至56%、55%,可以看出降雨达到5 h时,低危险区域面积占比下降幅度较大,高危险区域面积占比显著增加并随着降雨时间增加而逐渐增大。中危险区域面积占比出现了先增加后降低1%的情况,归结原因是降雨5 h后到t=10 h时仅有极少部分的滑坡低危险程度区域转为了滑坡中危险程度区域,而中危险程度区域转为高危险程度区域的面积相对较多。

本文利用既定的研究区参数(黏聚力、内摩擦角、土体重度),分析了坡度、地下水位深度(入渗深度)及PGA三种因素对滑坡位移值的影响(图9)。

图9(a)中为无降雨情况下,不同坡度(15°、20°、30°、35°)对应计算出的累积位移曲线图,在坡度值为15°、PGA为0.1g时所计算的位移值为负数,表明坡体不会发生滑坡,本文中取位移值为-0.01,即坡体稳定。由图9(a)可知:(1)在PGA值为0.2g时坡度值15°对应的位移值最小(0.44 cm),坡度值35°对应的位移值最大(1.63 cm);(2)PGA为0.2g时,不同坡度所对应的累计位移值相差较大,而PGA为0.5g时不同坡度所对应的累计位移值差距较小,主要原因是当坡度值较小时,PGA值越大对累计位移的影响度越大,坡度影响度随之减小。

图9 影响因子对滑坡位移影响图Fig.9 Influence diagram of influencing factors on landslide displacement

图9(b)为降雨情况下,降雨入渗深度Hw=0.6 m时不同坡度下计算出的累计位移,结果表明:(1) 降雨情况下,PGA=0.2g时,坡度为15°时对应的累计位移同样为最小值(1.31 cm),坡度值35°对应的位移值同样为最大值(2.12 cm);(2)在坡度为15°、PGA分别为0.2g、0.5g时的累计位移相比于图9(a)中无降雨情况下分别增加0.87 cm、0.84 cm,增加幅度分别为198%、30%,表明当坡度与PGA等其他参数值相同时,降雨入渗作用将增加滑坡位移值,当PGA值与坡度值较小时,增加幅度大;当PGA值与坡度值较大时,增加幅度小。

图9(c)列出了当坡度值30°,不同入渗深度(0.2 m、0.4 m、0.6 m、0.8 m)下的累计位移值曲线。结果表明:(1)PGA为0.1g时,入渗深度为0.2 m所对应的累计位移值最小(1.65 cm),入渗深度为0.8 m所对应的位移值最大(1.82 cm),增加幅度为10.3%;(2)当PGA值为0.5g时,四种入渗深度所计算出的累积位移值几乎相同。主要原因是坡度与PGA达到较高数值时,降雨入渗因素对累积位移的影响程度较小。

4 结论

本文结合GA降雨模型与Newmark累积位移模型在降雨地震耦合作用下对安全系数Fs进行了推导,对鲁甸县某区域开展无降雨情况、降雨无积水与降雨积水三种情况地震滑坡危险性预测分析,并对滑坡位移值进行坡度、降雨入渗及PGA的因素分析。

(1) 相对无降雨情况,加入降雨因素后,随着降雨时间的增加,滑坡低危险区域逐渐减少转变为滑坡中危险区域,地震滑坡高危险程度区域增加,计算区域中滑坡高危险程度区域占比上升幅度为1 100%,滑坡低危险程度区域占比下降幅度为35%,说明降雨与地震耦合作用对诱发滑坡的影响不容忽视,降雨较大的增加了地震滑坡的可能性。

(2) 当坡度较小时,PGA值对滑坡累积位移影响相比坡度影响较大,当坡度增加至一定数值时,PGA对滑坡位移影响程度远小于坡度;与无降雨情况相比,降雨入渗将使滑坡位移值增大,入渗深度在PGA值较小时对滑坡位移值影响较大,而在PGA值较大时对滑坡位移值相对影响较小。

(3) 降雨与地震耦合作用下的计算模型与降雨强度q、降雨时间t、岩石物理力学参数以及Dem数据的精度有较大关系,在进行耦合作用下滑坡预测分析时,获取精确的计算参数将极大提高该方法预测的精确度。