0.   引言

Newmark累积位移模型在地震滑坡风险评估中得到广泛应用，且不断得到方法上的改进和完善。Newmark建立刚塑性体的滑块模型，依据地震动峰值加速度值和地震力来计算滑坡位移，这也为基于GIS对地震滑坡的危险性评价提供了理论基础^[1]。基于典型模型方法，Wilson等对美国加利福利亚州(State of California)的洛杉矶地区(Los Angeles District)自然斜坡开展了区域地震滑坡的风险评估，对此方法的推广应用很有启发性^[2~3]；Wieczorek等对加利福利亚周圣马特奥镇(SanMateo)进行了地震滑坡易发性的实验性编图研究^[4]；Jibson等对加利福利亚州北岭市(Northridge)1994年地震诱发的滑坡进行了危险性评估^[5]。Roberto和Refice等利用简化Newmark模型与蒙特卡洛(Monte Carlo)方法相结合，在意大利地区进行了地震滑坡危险性评估^[6~7]；Godt等对汶川地震滑坡易发性开展了快速评估的尝试^[8]，结果表明此方法适用于较大区域(小比例尺)地震滑坡的评估。

国内的张春山等对夏呀河滑坡开展了单体滑坡的风险评估，推动了对单体滑坡灾害风险评估的方法和定量评估模式^[9]；李典庆等开展了考虑时间效应的风险评估，着重从时间效应风险滑坡风险评估的因素^[10]；乔建平等对区域滑坡进行了风险性概率分区^[11]；王涛等以M_s8.0级汶川地震为例，对震中附近区域使用简化Newmark位移模型法进行了区域性地震滑坡危险性快速评估，显示了评估方法的可行性^[12]；徐光兴等建立了地震强震动下边坡累积位移预测模型^[13]；陈晓利等则对庐山地震使用Newmark法进行了地震滑坡的预测研究，表明此方法具有可行性和实用性^[14]。

现阶段的Newmark累积位移法模型在滑坡危险性评估中存在以下不足：(1) 地震波加速度时程分为UP、NS、EW三个方向，使用回归方程取代单一方向相应的积分运算，考虑了震源或断裂的方向性，但是忽略了研究区斜坡的方向性；(2) 栅格计算的时候以单元栅格独立计算，忽略了相邻单元之间的关系，栅格与栅格之间未能搭建起关系，计算的结果存在偏差；(3) 在累积位移的计算中，无法做到单独栅格对某一地震加速度数据开展直接积分运算，降低了计算结果的准确性。本文尝试在野外调查数据的基础上，通过地理空间数据的统计分析，首次引入斜坡坡向和震源方位角的关系，对天水市北山开展基于Newmark累积位移模型的危险性分析，进一步完善地震滑坡的评估方法，为天水市减灾防灾提供理论基础和技术指导。

1.   天水市北山概况

天水市北山(下文简称“北山”)位于甘肃省天水市秦州区和麦积区交界，地处陇山山地与陇西黄土高原的过渡地带，是典型的黄土梁峁地貌，南北边界分别为藉河和渭河河谷冲积平原，藉河左岸冲积平原比渭河宽。北山东部为黄土塬，中部和西部为黄土梁，整体而言北坡比南坡坡度大。藉河在北山东南的社棠镇西侧汇入渭河。北山海拔最高点(1500 m)位于东部的黄土塬上，最低点(1120 m)位于渭河最东段河床上，高差380 m，北山总面积约41 km²。出露地层主要为第四系黄土，约占86%，厚度约15~45 m，南坡水眼寨附近出露面积约10%的新近系泥岩，北坡渭河岸坡出露少量下古生界牛头河群变质岩。北山主要的斜坡类型有：岩土二元结构斜坡、土质斜坡和岩土混合斜坡三类。天水市成纪新城选址于南坡藉河左岸河谷盆地(见图 1)。

图  1  天水市北山位置和滑坡分布图

Figure  1.  Location and landslide distribution of Beishan, Tianshui

下载: 全尺寸图片 幻灯片

北山位于青藏高原东缘，北部、东北部发育六盘山断裂和海原断裂，走向为北偏西。西南部为迭部—白龙江断裂，基本走向为北西西。南部有北东东向的成县盆地北缘断裂、平武—青川断裂和近东西走向的成县盆地南缘断裂。东部有走向近东西的渭河断裂和秦岭北麓断裂。北山区域属于潜在强震区或强震扰动区，周边区域历史上共发生55次M_s=6.0及以上地震，其中发生在麦积区范围的地震有734年天水7.0级地震和1654年天水南8.0级地震，震中位置分别为马跑泉镇一带^[15]和天水镇西北一带^[16](见图 2)。地震在该地区诱发了大量的黄土滑坡，据统计仅北山就发育了55个地震黄土滑坡(见图 1)，其中多为大型—巨型滑坡远程滑坡，平均滑距达969 m，地震滑坡面积约占北山面积的17.3%，给当地的生命财产安全造成了严重威胁^[17]。

图  2  天水市及周边断裂和历史地震(M_s≥6.0级)分布图

F1—西秦岭北缘断裂；F2—迭部—白龙江断裂；F3—礼县—罗家堡断裂；F4—六盘山断裂；F5—海原断裂

Figure  2.  Distribution of faults and historical earthquake sites (M_s≥6.0) in Tianshui City and its vicinity

下载: 全尺寸图片 幻灯片

2.   引入坡向角的Newmark累积位移法计算原理

Newmark模型基于极限平衡理论，该模型中滑块的永久位移是在地震加速度荷载作用下，斜坡块体沿着最危险滑动面发生瞬时失稳后位移不断累积所得。当施加于最危险滑动面处的加速度超过临界加速度时，块体即沿滑面发生滑动，将外荷载加速度与临界加速度的差值部分曲线对时间进行二次积分即可得到永久位移值。

临界加速度a_c指在地震动荷载作用的条件下，斜坡滑体的下滑力等于斜坡内的抗滑力时(即力的极限平衡状态)的地震动加速度值在斜坡坡向角的分解值a′_t(见图 3)。由于一次地震波的方向是固定的，而斜坡的坡向值是不固定的，无法对栅格进行定量计算，通过数学变换的方法，斜坡临界加速度分解到单一固定地震波方向上，此时的斜坡临界加速度为a′_c，计算结果与地震波分解到斜坡方向上是一致的。通过比较静力和地震动力条件下滑体的受力状态，利用无限斜坡法计算安全系数Fs间接公式(1) 推导得出，g为地球重力加速度，α为斜坡面的倾角,β为斜坡坡向角,θ为地震波加速度方位角：

图  3  Newmark模型滑体静力受力图和地震加速度水平分解图

Figure  3.  Static force diagram of Newmark model and diagram of seismic acceleration horizontal decomposition

下载: 全尺寸图片 幻灯片

其中Fs为无限斜坡法计算出来的安全系数，其计算如式(2) 所示：φ′为斜坡土体的内摩擦角，α为斜坡的坡角，c′为斜坡的内凝聚力，γ为坡体物质密度，γ_w为水的密度，m为滑体条块被水浸入的厚度比例，t为坡体的厚度。

利用临界加速度a′_c计算滑块累积位移，还需要地震动加速时程曲线数据，加速度值小于等于a′_c的部分对滑块不产生位移，而将大于a′_c的加速度时程对t进行一次积分求出滑块对应时间的速度时程曲线见图 4，按式(3) 二次积分可得滑体累积位移时程曲线：

图  4  考虑斜坡坡向角的Newmark累积位移求解原理图

Figure  4.  Principle diagram of Newmark cumulative displacement model in consideration of the slope angle

下载: 全尺寸图片 幻灯片

3.   北山地震黄土滑坡易发性评价

根据地层岩性和岩土体工程特性，按照工程地质分类标准^[18]将研究区内工程岩组划分为4种类型：坚硬岩组、较坚硬岩组、较软弱岩组、软弱岩组。较软弱岩组和软弱岩组则均为第四纪黄土松散沉积物，其土体参数值基于北山的红旗山滑坡和渭河北岸的锻压机床厂滑坡两处岩土体测试结果，再综合该地区的已有成果资料，得到本文研究各类岩土体的物理力学参数(见表 1)。

表  1  北山岩土体物理力学参数取值表

Table  1.  Physical and mechanical parameters of rock and soil in Beishan

岩组	坚硬 (变质岩)	较坚硬 (砂砾岩)	较软(松散 沉积物， 藉河， 红旗山)	软弱(松散 堆积物，距河 流较近， 锻压厂)
φ′/°	40	35	24.5	22.2
c′/Pa	220000	140000	90400	59500
γ/(N·m -3)	24500	23000	17800	17500

下载: 导出CSV 

| 显示表格

利用DEM数据生成斜坡坡度栅格图，结合研究区滑坡滑动面倾角进行插值改进，得到斜坡潜在滑动面倾角(α)栅格图层(见图 5a)；根据1:5万地质图基础数据，进行数据矢量化，然后利用图幅内的红旗山滑坡和渭河北岸的锻压机床厂滑坡的岩土物理力学测试结果进行属性赋值，最后转为有效内摩擦角(φ')、内聚力(c)和岩土体重度(γ)栅格图层(见图 5b—图 5d)；北山共发育108个新老滑坡，通过详细野外调查，获取了每一个滑坡的滑坡体大致厚度数据，利用自然邻域法插值分析，从而得到北山土体的厚度(t)分布栅格数据(见图 5e)；收集麦积区内6个村雨量站在2013年6~7月两次强降雨中记录的过程雨量，取两个月的平均值为点源数据，利用ArcGis进行降雨量插值，生成过程雨量均一化数据后作为土体饱和系数(m)的栅格数据(见图 5f)；使用DEM提取栅格的坡向角进行分析，得到优势坡向角值θ (见图 5g)。

图  5  Fs求解参数因子栅格图

Figure  5.  Grid diagramaboutparameters factor of Fs

下载: 全尺寸图片 幻灯片

由公式(2) 可以求出北山区域静态安全系数Fs的栅格图层，将计算结果分为六个等级(第六级为坡度小于5°的区域)，得到不同分级的静态安全系数Fs的分区图(见图 6)。研究区东侧的水眼寨、叶家沽沱、董家咀、孙家坪以及红旗山地区的安全系数最低，约占研究区面积的4%；周家山、小见河、罗圈湾、东胡家沟到孙家坪等周边的Fs为1~2，属于较敏感区域，占研究区面积的23.5%；南部的藉河河谷盆地、东侧的黄土塬和北部的渭河河谷较为平缓，坡角在5°范围，安全系数极高，面积约为38%；其他地方安全系数处于2~5之间，约占19.3%。

图  6  北山静态安全系数Fs栅格图

Figure  6.  Grid diagram of static safety factor Fs in Beishan

下载: 全尺寸图片 幻灯片

利用公式(1) 计算得到优势坡向角下临界加速度栅格图(见图 7)和劣势坡向角下的临界加速度栅格图(见图 8)，其中优势坡向角指的是地震波和斜坡坡向角方向一致，劣势坡向角指的是地震波与斜坡坡向角垂直。临界加速度值分为五个等级：极高敏感度(0~0.2g)，高敏感(0.2~0.4g)，中敏感(0.4~0.6g)，低敏感(0.6~1.0g)和极低敏感(>1.0 g)。对比图 7和图 8可知：优势坡向角下极高敏感度明显比劣势坡向角下的范围大，二者的面积比例分别为37.8%和10.6%；高敏感度面积比分别为17.8%和11.4%。

图  7  优势坡向角下的斜坡临界加速度a′_c1

Figure  7.  Critical acceleration a′_c1 of slope under dominant slope angle

下载: 全尺寸图片 幻灯片

图  8  劣势坡向角下的斜坡临界加速度a′_c2

Figure  8.  Critical acceleration a′_c2 of slope under adverse slope angle

下载: 全尺寸图片 幻灯片

4.   北山地震黄土滑坡危险性分析

徐光兴等利用汶川地震时记录到的大量强震动记录，建立了基于临界加速度比、Arias强度和地震剩余强度的永久位移预测模型^[13](式4)。根据根据第五代《中国地震参数区划图》(GB 18306-2015)，北山区域PGA取0.3g，计算得到优势坡向角下的累积位移Dn₁栅格图和劣势坡向角下的累积位移Dn₂栅格图(见图 9)，分为五个级别：极高危险、高危险、中等危险、低危险和极低危险。

图  9  优势坡向角和劣势坡向角的斜坡累积位移栅格图

(a)优势坡向角(b)劣势坡向角

Figure  9.  Grid diagram of accumulative displacement under dominant slope angle & adverse slope angle

下载: 全尺寸图片 幻灯片

Newmark模型求出的累积位移Dn与边坡失稳不存在直接联系，但可将其与实际的地震黄土滑坡编目进行对比分析以确定Dn与地震黄土斜坡失稳概率P(Dn)之间的关系。据Jibson与葛华的相关研究表明^{[15, 19]}，以地震诱发滑坡编目为数据源，建立的地震黄土斜坡失稳概率指数K与累积位移Dn的回归方程如公式(5)：

其中中括号内为Weibull方程，P(D_n)表示斜坡在地震动加速度作用下产生累积位移Dn时诱发滑坡的失稳概率，系数k表示P(Dn)最大值，而a与b作为选定参数，用每个分组内处于滑坡区域的累积位移Dn栅格单元数量除以分组内Dn栅格数量的值算出。使用ArcGIS进行重分类(Dn取中间值)，计算出每一类位移值的总面积A，再将地震黄土滑坡数据与累积位移数据重分类图层进行相互分析，计算出累积位移内的地震黄土滑坡面积a，取P(Dn)=a/A，见表 2。

表  2  累积位移Dn的Weibull函数拟合参数表

Table  2.  The table of Weibull function fitting parameter of cumulative displacement

Dn	P(Dn)
0.005	0.0891
0.05	0.0936
0.5	0.1566
1.5	0.1720
2.0	0.1921

下载: 导出CSV 

| 显示表格

通过公式(5) 和表 2的参数，对累积位移栅格数据进行计算，得到优势坡向角和劣势坡向角下的斜坡发生失稳概率分布图。使用自然断点分级法(Jenks)把发生斜坡失稳概率分为五级：极高、高、中等、低和极低。图 10为优势坡向角下的斜坡失稳概率栅格图，最大失稳概率为19.2%，最小为0.1%，与累积位移分区图(见图 9)基本吻合：水眼寨、叶家沽沱、陈家山、周家山、红旗山、缑家窑和董家塬一带失稳概率最高；藉河和渭河河谷、河漫滩及河谷冲积平原地带失稳概率最小，在2%以下；县家路村—赵集寨—坡阳里一带的斜坡失稳概率较低。

图  10  优势坡向角下斜坡失稳概率图

Figure  10.  Diagram of instability probability of slope under dominant slope angle

下载: 全尺寸图片 幻灯片

图 11为劣势坡向角下斜坡失稳概率图，优势坡向角下的失稳概率区域面积明显比劣势坡向角下失稳概率区域面积大。对比二者的失稳概率图可以看到：失稳概率极高和高的区域明显变小，由32.4%下降到10.1%；周家山—孙家坪北、水眼寨、叶家沽沱和陈家山—西胡家湾一带地区的失稳概率最大值为17.6%，失稳概率为高—极高。根据北山现有滑坡分布和不稳定斜坡发育的调查结果显示：优势坡向角下的计算结果和北山不稳定斜坡调查结果基本吻合。成纪新城西区闫家河村到孙家坪一带和东区水眼寨村南一带为失稳概率高区域，存在较大地震黄土滑坡风险，建议加强该地区滑坡变形监测，开展减灾应对措施。

图  11  劣势坡向角下斜坡失稳概率图

Figure  11.  Diagram of instability probability of slope under adverse slope angle

下载: 全尺寸图片 幻灯片

5.   结论

本文通过天水市北山野外调查和室内分析相结合的方法，开展基于Newmark累积位移法的地震黄土滑坡危险性分析。分析了研究区的孕灾背景，结合现有滑坡和不稳定斜坡的调查数据，在充分考虑了震源方位角和斜坡坡向角的角度关系基础上，对优势坡向角和劣势坡向角下的斜坡易发性进行计算和对比，然后用回归方程求出PGA为0.3g时斜坡的累积位移，最后计算斜坡在地震荷载下的失稳概率，得到结论以下：

(1) 天水市北山及周边构造发育、历史地震频发，诱发了大量地震黄土滑坡，滑坡面积占到北山面积的27%，该区域属于潜在强震区或强震扰动区。

(2) Newmark累积位移模型中引入研究区优势坡向角变量，对地震黄土滑坡危险性评价的结果准确性有提高作用。通过北山的危险性分析结果对比，在不考虑斜坡优势方位角下北山的斜坡失稳概率分区图的识别率低，且与调查结果不符；而劣势坡向角下的结果和实际情况存在较大差距；优势坡向角的斜坡失稳概率和该区域调查结果吻合，在1:5000比例尺区域有较好实用性和合理性。对Newmark模型在地震黄土滑坡危险性评估中引入斜坡坡向角和震源方位角关系，有利于提高小区域(大比例尺)风险评估的准确性。

(3) 北山的周家山—孙家坪一带、水眼寨一带、叶家沽沱一带和罗圈湾—陈家山一带的在地震扰动下斜坡失稳概率为16.9%~19.2%，为极高危险区。成纪新城西区的红旗山—孙家坪一带、东区的水眼寨南一带均为高—极高危险区，建议加强对该地区地质灾害监测，开展减灾设施建设，保证在建和拟建工程的安全运营。