考虑基质吸力作用的Newmark改进模型在地震滑坡风险评价中的应用

冯 卫，唐亚明，赵法锁，陈新建

(1.长安大学地质工程与测绘学院，陕西 西安 710064；2.中国地质调查局西安地质调查中心，陕西 西安 710054)

对于地震滑坡而言，地震滑坡危险性分析和风险区划是降低财产损失和保护生命的有效手段之一，并逐渐成为研究热点之一。Newmark[1]于1965年提出了在地震力作用下通过累积位移来分析边坡变形的方法，对在地震作用下边坡的稳定性进行了定量研究，该模型对地震诱发的浅层刚性块体发生平移或旋转滑动的永久位移评估适用性较好[2-3]。众多国外学者相继利用Newmark模型开展地震滑坡的危险性评价和风险评估[4-6]，编制了一系列危险性评估图件。国内利用Newmark模型进行区域地震滑坡危险性研究工作起步较晚，王涛等[7]以汶川Ms8.0级地震重灾区的11县市为例，提出了基于简化Newmark位移模型的地震滑坡危险性应急快速评估方法。近几年来，相继有部分学者对该模型进行改进[8-11]，不断提高了模型的准确性和实用性。

值得注意的是，现阶段的Newmark评估模型主要针对饱和土的边坡稳定性分析，对于地下水位以上由负孔隙水压力提供的部分抗剪强度则忽略不计，这对于滑动面的主要部分处在地下水位以下时较为合理，但对地下水位很深或考虑出现浅层滑动的坡体，基质吸力作用不容忽视[12]。大量研究表明[13-19]，非饱和土基质吸力的有效应力作用，使得边坡的安全系数增大，这对边坡稳定性分析来说具有重要的现实意义。本文将Newmark模型扩展至非饱和土力学领域，考虑了非饱和土的基质吸力作用，对静态安全系数公式进行改进，进一步完善了地震滑坡评估方法，并针对甘肃礼县幅区域开展地震滑坡危险性分析和风险评价。

1 考虑基质吸力作用的Newmark累积位移模型

1.1 传统Newmark累积位移模型

Newmark模型[1]基于极限平衡理论，将滑体视为刚体，假设滑体内部不发生形变，当受到的外力作用小于临界加速度时，坡体不产生永久位移；当受到的外力作用大于临界加速度时，则会产生有限位移。该模型通过计算滑体在地震动加速度作用过程中累积的永久位移来评价斜坡的稳定性。

目前，在计算斜坡稳定性时，普遍采用考虑降雨入渗影响的安全系数公式，其受力状态见图1。

(1)

式中：φ′——有效内摩擦角/(°)；

c′——有效黏聚力/kPa；

α——坡角/(°)；

γ——坡体物质重度/(kN·m-3)；

γw——水的重度/(kN·m-3)；

t——坡体厚度/ m；

m——破坏面以上的地下水埋深与破坏面深度的比值。

图1 Newmark模型静态受力分析示意图Fig.1 Static force analysis diagram of Newmark model

临界加速度表示为：

ac=(Fs-1)gsinα

(2)

式中：g——重力加速度/(m·s-2)。

Newmark位移计算采用Jibson和Harp等研究得到的经验公式：

lgDn=1.521lgIa-1.993lgac-1.546

(3)

式中：Dn——Newmark位移/ m；

Ia——地震强度数据/(m·s-1)。

Ia值的通用公式为：

(4)

式中：a(t)——强震记录中单一分量的地震动加速度时程序列/(m·s-2)；

Td——地震动加速度记录的总持续时间/s；

t——时间/s。

1.2 改进后的Newmark累积位移模型

笔者在总结现有模型的基础上，认为现阶段Newmark评估模型中的安全系数在计算上存在一定的局限性，该公式未考虑地下水位以上由负孔隙水压力提供的部分抗剪强度，这对于滑动面的主要部分处在地下水位以下时较为合理，但对地下水位很深或考虑出现浅层滑动的坡体，基质吸力作用不容忽视。而地震形成的滑坡中，存在大量浅表层滑坡，且滑坡大部分位于地下水位以下，因此有必要对该安全系数公式进行适当改进。

根据Fredlund提出的非饱和土抗剪强度公式可知：

τf=c′+(σ-ua)tanφ′+(ua-uw)tanφb

(5)

式中：σ-ua——破坏时破坏面上的净法向应力/kPa；

ua-uw——基质吸力/ kPa；

φb——抗剪强度随基质吸力增加的速率/(°)。

本文结合非饱和土的抗剪强度公式，将基质吸力作用考虑在内，对Newmark模型的静态安全系数公式进行了改进，将其表示为：

(6)

式(6)中的参数在前文中均已出现且含义完全相同，此处不再重复说明。

2 算例分析

2.1 研究区环境地质概况

选取1∶5万标准图幅中的甘肃礼县幅(I48E011013)为研究区(图2)。该研究区位于甘肃省陇南市礼县北部，东侧区域地层岩性以新近系红黏土、第四系中晚更新世风积黄土为主，西侧区域以燕山期中细粒花岗岩、泥盆系板岩、砂岩及石炭系炭质板岩、炭质页岩为主。构造上地处南北构造带与东西构造带的构造结内，地质构造复杂，地震活动性强。1654年甘肃罗家堡8级地震震中距研究区直线距离约30 km，1879年甘肃武都南8级地震震中距研究区直线距离约90 km。

图2 研究区地理位置图Fig.2 Geographical location map of the study area

2.2 参数选取

在选取工程地质岩组时，考虑到Newmark方法的实用性，对研究区的地层岩性进行了简化，将研究区的岩土体划分为5个岩组，即坚硬岩组(Ⅰ)、较硬岩组(Ⅱ)、较软岩组(Ⅲ)、软弱岩组(Ⅳ)和松散岩组(Ⅴ)(图3)。

图3 工程地质岩组图Fig.3 Map of engineering geological rock group

坚硬岩组主要是黑云母斑状花岗岩、中细粒花岗岩，结构较完整，斜坡稳定性好；较硬岩组主要是中泥盆统舒家坝组和西汉水组浅变质砂岩、灰绿色板岩、粉砂质板岩、粉砂岩，结构较完整，斜坡稳定性较好；较软岩组主要是中石炭统灰黑色炭质板岩、泥质灰岩、炭质粉砂岩、炭质页岩，岩体强度较低，斜坡稳定性一般；软弱岩组主要是新近系红黏土，该岩组成岩性差，较松软，呈披覆状不整合于基底之上，可与上部黄土构成相对隔水层而引发滑坡；松散岩组以中晚更新世风积黄土为主，披覆于新近系红黏土之上，是该地区土体滑坡的主要发育层位。因松散岩组和软弱岩组中存在非饱和区，故需考虑基质吸力作用。本文基于区域性钻孔取样和探槽取样测试结果，并综合该地区已有的成果资料，最终得到研究区所用的各类工程地质岩组力学参数(表1)及松散岩组和软弱岩组的土水特征曲线(图4)。

表1 工程地质岩组力学参数取值

图4 土水特征曲线Fig.4 Soil water characteristic curve

在选取地形坡度时，利用现有的1∶5万数字高程模型(DEM数据)，通过计算相邻栅格单元之间高差与水平距离之比获得。野外调查发现，地形坡度小于10°的斜坡基本处于稳定状态，发生较大规模滑坡的可能性很小，为进一步提升评估计算效率，对相应区域斜坡不予计算(图5)。

图5 地形坡度图Fig.5 Map of topographic slope

2.3 易发性分区

分析滑坡稳定性时，滑体状态按照干燥状态考虑，取破坏面以上的地下水埋深与破坏面深度的比值(m)为0，由于研究区大多数滑坡破坏面深度较浅，这里滑体厚度(t)近似取5 m，分别利用传统模型和改进模型，计算区域斜坡体静态安全系数Fs分布，进而求出相应的区域临界加速度分布。由于临界加速度表征了因斜坡固有属性导致地震诱发坡体失稳的潜势，可以作为区域地震滑坡易发性评估的依据[7]。本文采用突变点法，经过统计分析，从中找出突变点作为易发程度分区界线值，将临界加速度分为极高易发区(小于0.16g)、高易发区[0.16g，0.46g)、中易发区[0.46g，0.72g)、低易发区[0.72g，1.04g]和极低易发区(大于1.04g)五个等级(图6)。

图6 易发性分区图Fig.6 Division map of susceptibility

根据模型改进公式(6)，由于松散岩组和软弱岩组中考虑了基质吸力作用，使得这部分区域内计算的静态安全系数值比改进前显著提高，相应计算出的区域临界加速度值也相应增大，进而导致易发程度相应降低。由图6可以明显看出这部分易发性分区发生了相应变化。而其他岩组区域由于计算参数未发生改变，故模型改进前后易发性分区无变化。

2.4 危险性分区

危险性分区是在计算区域静态安全系数和区域临界加速度的基础上，开展Arias强度计算。由于区内暂未搜集到地震动加速度时程序列等相关参数，本文计算Arias强度利用Roberto R分析意大利17次强震的190条加速度时程记录拟合的经验公式[20]：

Ia=0.004(PGA)1.668

(7)

式中：PGA——地震动峰值加速度/(m·s-2)。

根据中国地震局兰州地震研究所编制的对应50年超越概率10%的地震动参数区划，求得相应的Ia分布情况，进而求出模型改进前后的50年超越概率10%条件下的Newmark累积位移值，依据Newmark累积位移值的大小将地震滑坡的危险性等级划分为高危险区(大于1 cm)、较高危险区(0.1，1 cm]、中危险区(0.05，0.1 cm]、较低危险区(0.02，0.05 cm]和低危险区(0，0.02 cm]五级(图7)。

由于危险性分区是在易发性分区基础上开展的，模型改进后同样导致松散岩组和软弱岩组区域内的危险性相应降低。从图7可以明显看出这部分危险性分区发生了相应变化，而其他岩组区域由于计算参数未发生改变，故模型改进前后危险性分区无变化。

图7 50年超越概率10%条件下的地震滑坡危险性分区图Fig.7 Hazard division map of earthquake landslides under 10% probability in 50 years

2.5 风险评价

风险评价是在Newmark法计算得出斜坡体累积位移分布结果(危险性结果)的基础上，对危害后果(人员、财产及易损性)进行量化(表2)，将危险性量化结果与危害后果量化结果(危害后果=人口数量或财产价值×易损率)进行栅格乘积运算，分别得出传统模型和改进模型50年超越概率10%条件下的地震滑坡风险评价结果(图8)。

风险评价是在危险性分区基础上开展，利用改进模型计算使得松散岩组和软弱岩组区域内的风险程度相应降低。从图8可以看出这部分风险分区情况发生了相应变化，而其他岩组区域由于计算参数未发生改变，故模型改进前后风险分区无变化。

表2 承灾体类型、价值及易损性估算值

图8 50年超越概率10%条件下的地震滑坡风险评价图Fig.8 Risk assessment map of earthquake landslide under 10% probability in 50 years

2.6 数据验证

为了分析和比较模型改进前后所得到的地震滑坡危险性和风险评价结果，本文采用受试者工作特征曲线(Receiver Operating Characteristic，ROC)对计算结果进行验证。根据ROC曲线下的面积值(简称AUC值，Area Under the Curve)作为衡量准确度的指标。一般认为，AUC值在0.7～0.9之间表示结果具有一定的可靠性，其值越接近1可靠性越高[21]。通过计算发现，改进模型所得易发区和危险区的AUC值明显高于传统模型计算得出的AUC值，说明改进效果显著。而改进模型所得风险区的AUC值虽然没有明显提高，但仍高于传统模型计算出的AUC值(图9和表3)。总之，改进模型的计算结果更为可靠合理，具有较高的实用价值。

表3 传统模型和改进模型AUC值

图9 传统模型和改进模型ROC曲线对比图Fig.9 Comparison diagram of ROC curve between traditional model and improved model

3 结论

(1)针对浅表层土质滑坡，结合非饱和土的抗剪强度公式，将基质吸力作用考虑在内，改进了Newmark模型的静态安全系数表达式。

(2)对甘肃礼县幅区域内的地震滑坡开展易发性分区、危险性分区及风险评价的过程中，利用Newmark改进模型计算所得结果在松散岩组和软弱岩组区域内均得到不同程度的优化，这主要由于松散岩组和软弱岩组中考虑了基质吸力作用，使得这部分区域内计算的静态安全系数值比改进前显著提高，进而使得易发程度、危险程度及风险程度相应降低。而其他岩组区域由于计算参数未发生改变，故模型改进前后相应分区无变化。

(3)提出的考虑基质吸力作用的Newmark累积位移模型，弥补了原公式中忽略基质吸力作用的缺陷，使得计算结果更为合理，研究成果可为相关地区开展地震滑坡易发性分区、危险性分析和风险评价提供新的参考。