0 引言

大量震害调查表明，地震诱发的液化灾害是典型的地震灾害之一.1964年日本新瀉地震、美国阿拉斯加地震中大面积的液化，使地基承载强度丧失，大量的建筑物倒塌，道路、桥梁以及地下管道等基础设施都会遭到严重破坏¹⁾，1976年中国唐山大地震、1999年中国台湾集集地震(Chu et al., 2004)、2008年汶川地震(曹振中等，2010)、2011年新西兰地震(Cubrinovski and Robinson, 2016)等都显示了液化灾害的巨大破坏性.

1) Jibson R W, Michael J A. 2009. Map showing seismic landslide hazards in Anchorage, Alaska. US Department of the Interior, US Geological Survey.

减轻土壤液化灾害的首要手段是能科学合理地预估和预判工程场地的液化情况，主要工作包括两部分，一是区域性的土壤地震液化风险估计与区划，二是对重要工程场地详勘及其液化判别，二者目标不同，不能互代，但可互相补充.对重要工程场地详勘及其液化判别，主要依赖于钻孔和现场测试，这方面我国抗震设计规范中已经有明确的规定.但是，这种强烈依赖于钻孔及其实测资料的方法，难以适应工作区尺度较大的情况.

液化风险评估和区划研究起源于19世纪70年代，一些先进国家已将较大尺度液化灾害风险评估和区划工作纳入了国家制图计划项目中，美国和日本等多个国家已经编制了有关液化的各种形式的区划图(Youd and Perkins, 1978, 1987).美国国家地震减灾计划(NEHRP)早在20世纪70—80年代就制定了绘制液化灾害区划图的计划，还为此专门成立各个州的液化咨询小组，内容包括了大范围地区液化灾害评估工作，如基于地质准则绘制了旧金山9个区面积为1000多km²的比例尺1:100万大尺度的液化敏感度图，给出了这些地区的液化风险估计.2016年发生的中国台湾高雄地震，土壤液化严重，引起中国台湾民众强烈关注.随后中国台湾“行政院”立即组织研发液化潜势评估查询系统，历经一个月的时间发布了台北市、新北市地区1:15万液化风险图，民众响应热烈(新浪网，2016).

一些学者根据地质和地貌条件研究了液化的触发条件(李程程，2013).Knudsen和Bott研究发现，液化场地通常在地质年代、沉积环境、坡度、海拔以及到水体的距离等具有较为统一的属性²⁾.通过对历史液化数据的研究，也发现了地质年代、沉积环境、水文条件等与液化发生很强的关联性.NEHRP授权绘制的1:100万旧金山湾地区液化区划图，是依据地质准则和较为详尽工程地质图并依照经验判断完成的.其中，最重要资料来源是美国拥有的1:24000的地质图以及旧金山湾地区详尽的工程地质资料，但这两点目前在我们国家还难以做到.同时，对网格中各节点是否液化的判断，没有提出预测模型，而是从以往经验得到表格方式中获取，而且对此方法也没有验证.

2) Knudsen K L, DeLisle M D, Clahan K B, et al. 2000. Applicability of Quaternary geologic mapping in assessing earthquake-induced liquefaction hazard: San Francisco Bay Region: Sixth International Conference on Seismic Zonation, Earthquake Engineering Research Institute.

对于区域性土壤地震液化风险估计与区划技术，我国大陆研究水平还存在较大差距，远远不能适应我国经济建设发展需求.我国国土辽阔，地质环境介质复杂，三大城市群和众多大城市地震形势严峻且位于可液化场地之上.近期党中央国务院提出了加强我国自然灾害风险普查以及加强我国土壤液化风险防范能力的要求，发展适于我国情况的区域土壤液化风险评估方法就显得尤为迫切.另外，拥有土壤液化风险评估简化方法，就可以在地震后，通过实时发布的地震动参数，快速了解震区土壤液化可能的分布，对地震快速响应、损失估计和救援计划制定都很有帮助，但我国这两方面研究成果均为空白.

鉴于此，本文将总结以往成果，遴选先进的指标体系，以我国典型地震液化资料为背景，构建区域土壤液化预测模型，并采用近期大地震液化调查资料予以检验，以期为我国较大尺度区域的土壤液化预测提供实用方法，并为震后快速评估区域土壤液化震害提供直接的技术支持.

1 参数选取和数据样本

土壤液化是一个复杂的动力响应问题，影响因素较多，大致分为地震动大小和埋藏条件两类，并且可以采用一些一级代表参数用来预测给定地震动下土体液化的大致情况.

1.1 参数选取

首先，砂土密实程度是决定液化是否发生的首要因素，本文选用剪切波速V_s来描述土体密实程度.正如小应变下土体剪切模量G_max=ρV_s²所表达的那样，一般土体的密度ρ相差不大，则选用剪切波速V_s可直接用来估量土体的初始剪切模量.而研究表明(Seed et al., 1970)，砂土初始剪切模量与其相对密度密切相关.以往震害资料调查也表明，可液化土层的场地一般较为软弱，V_s一般较小，无论液化层多深，目前还没有发现实测剪切波速超过220 m·s^-1的非液化土层.同时，V_s作为一种有发展前途的场地性能衡量指标，与如标准贯入试验(SPT)和静力触探试验(CPT)有一定关联性，目前实验室内也可以进行土单元和土体的V_s测试并进行相关液化研究，因此V_s可以作为一个独立的参数表征土体的抗液化特性.同样重要的是，我国规范中已经有明文规定，工程勘察中要进行土层的剪切波速测试.而且，Wald(1996)、Wald和Allen(2007)统计得到了V_s30(地下30 m内土层平均剪切波速)估算方法，给出两套V_s30的经验关系，一是地形地貌显著的构造活跃区，二是地形变化缓和的大陆稳定区，进行得到了全球场地V_s30估算模型，并将之应用于全球地震快速评价响应系统(PAGER).研究表明，土体的剪切模量与V_s30直接相关，总体来说受孔隙比和平均有效应力的控制.但又由于平均有效应力对地下30 m的平均剪切波速的影响不太显著，因此我们认为孔隙比是影响剪切模量的重要因素，孔隙比又与土体的结构和粒径分布息息相关.剪切波速随着土体颗粒的增大而增大，坡度越低，沉积的颗粒越小.这样，V_s30就成为区域土壤液化预测模型中的第一个参数.这与本文致力于研究基于DEM的大尺度液化区划技术的快速评估的目标不谋而合，V_s30可以在USGS网站(http://earthquake.usgs.gov/hazards/apps/vs30/custom.php)网址上进行下载，为我们的研究提供了基础数据.

其次，土壤的饱水状态是决定是否发生的另一重要因素.土壤的饱水状态一直是区域土壤液化预测要素提取中的难题，以往只能通过详尽的水文地质勘察结果获取.NEHRP绘制旧金山湾地区液化区划图中，土壤的饱水状态就是通过详尽的水文地质勘察资料给出的，完成这一工作颇费周折.但是，因资料限制，这一方式难以推广到其他地区，土壤饱水状态的描述需要另寻蹊径.令人高兴的是，通过资料研究我们找到Beven和Kirkby(1979)提出的一种可以快速，方便和经济地体现水文特征的模型，它是以变源产流为基础的半分布式水文模型，即TOPMODEL(Topography-based hydrological model).该模型可以利用复合地形指数CTI，也称为湿度指数，可反映下垫面空间变化和对流域水文模型循环过程的影响，其表达为CTI=ln(a/tanβ)，其中a是通过每单位等高线长度的上游集水面积，β是相应的坡度.该参数将地貌特征与水文地质有机结合，具有明确的物理意义，并且验证结果良好.因此，本文采用这一先进的水文参数，作为土壤饱水状态指标.对于地形指数，首先要计算每单位等高线长度的上游集水面积a，然后计算流域中已知坡度的各个点的ln(a/tanβ)，可以将流域按最大坡度和等高线的交叉格网分割成大量的子域进行计算，以Lanshaw区域为例，计算结果与Kirkby(1976)制作的Lanshaw区域的饱和度分布图高度吻合.当CTI>10.0时，饱和度Sr>75%；9.0 < CTI≤10.0时，饱和度50% < Sr≤75%；8.0 < CTI≤9.0时，饱和度25% < Sr≤50%；CTI≤8.0时，0 < Sr≤25%.因此我们可以看到，地形指数CTI=ln(a/tanβ)的值与土体饱和度关系密切，可以作为饱和度的表征.本文选用分辨率为1 km的HYDRO1K DEM数据为水文分析的基础资料，HYDRO1K是由美国USGS的EROS Data Center地球资料观测系统数据中心与联合国环境计划——全球资源信息中心(UNEP/GRID, United Nations Environmental Program/Global Resources Information Database)联合制作的全球基于DEM的水文地形高程模型.它是由USGS发布的30arc-second(900 m)的分辨率的GTOPO30 DEM数据得到的，HYDRO1K包含一系列有关水文的矢量和栅格数据集，在GTOPO30数据基础上进行了河网校正以及洼地填充，提供了一套标准的地理参照数据集，对管理评估或处理水文信息的用户很有价值，可以免费在Internet上下载，覆盖全球(除格陵兰岛及极地以外)的所有陆地.

再者，土壤沉积年代与液化可能性密切相关.但直接获取沉积年代参数有一定困难，由于河流沉积物的搬运和沉积过程，离河流越远沉积层沉积年代就越久远.并且，到较大河流的距离DR也可以表示饱和度，因为到河流的距离越近地下水位越浅.因此，本文也采用到河流的距离DR作为区域土壤液化预测模型中第三个参数.

最后，土壤是否液化及液化规模与地震动强度密切相关.本文的区域土壤液化预测方法，目标是震前可估计各地区液化风险以及震后可快速评估震区的土壤液化情况，因此要求这一参数具备容易且快速获取的特点，本文采用地面峰值加速度(PGA)为参数.对于震前，我国新一代区划图已经对我国各地区提供了很详细的基本地震动的PGA.对于震后，本文可利用目前美国地质调查局Shake Map.Shake Map系统可对全球重大破坏性地震及时发布预测地震动图，可以在一个新的地震发生后，能够自动快速基于网络在一定比例的图上显示各个台站地震动情况，并且生成峰值地面加速度(PGA)等值图、峰值地面速度(PGV)等值图和峰值加速度反应谱等值图(0.3 s、1.0 s、3.0 s)等图形和数据产品，可用于我国震区土壤液化情况的快速评估.

1.2 数据样本

1976年7月28日唐山大地震液化举世闻名，地形地貌和水文地质特征具有典型性，并且通过国内外专家学者的不断研究，资料较为齐全(杜晨晓等，2010；刘晓霞等，2018)，对此次地震以及液化灾害等问题进行了详尽探讨与资料收集，为本文回归区域土壤液化预测模型参数打好了研究基础.

本文样本区域如图 1所示.该区域北依燕山，南濒渤海，地势北高南低.在地貌上，由于受纬向构造系的影响和控制而有很大的分异：北部为燕山及其余脉上升的低山丘陵区；中部是遭受轻微侵蚀的、准平原化的山前平原区；南部是广阔的滨海平原沉降区.

唐山地震后，中国建筑科学研究院对液化区进行了详细调查，利用航拍片判释结合现场测试对完成了震区液化综合评价，得到的液化区的分布如图 2所示.

本文将地震后三天的“7.28”航拍解译图运用ArcGIS进行栅格化，按100 m×100 m格网划分，样本区域如图 2所示，其中红色部分为液化样本，淡粉色为未液化样本.样本区域面积约9000 km²，共拥有896177个栅格点，包括218514个液化样本点，677663个非液化样本点，以此形成样本数据库.

对本文区域土液化预测模型中所选取的三类土壤特征以及地震动参数，包括剪切波速V_s30、综合地形指数CTI、场地到河流的距离DR和峰值加速度PGA.唐山样本数据库中四个参数的样本图形分别如图 3至图 6所示.样本区V_s30、CTI、DR和PGA的数据范围分布情况分别如图 7至图 10所示.

由图 3至图 10可见，样本中，对V_s30，非液化区范围为180~580 m·s^-1，液化区为180~345 m·s^-1；对CTI，非液化区多为0~10，液化区多为10~18；对DR，非液化场地为0~6.25 km，液化场地为0~4.25 km；对PGA，非液化区范围为0.35~1.1 g，液化区范围为0.35~0.75 g.从中可见，土壤几个主要特征参数的范围符合现有认识，如V_s30、DR越大液化发生可能性越小，CTI越大越容易发生液化.特别是数据库中，土壤几个主要特征参数的范围较宽，使得数据库具备较好的完备性.

2 计算公式及可靠性分析 2.1 区域土壤液化评估公式

运用GIS工具，将各参数代表值赋到每个点的属性值中.利用提取的896177个数据，采用二元Logistic回归，建立液化预测模型，得到区域土壤中场地液化评估表达式为

(1)

其中，Y为场地液化评估值，V_s30为场地平均剪切波速、CTI为综合地形指数、DR为到河流的距离以及PGA为峰值加速度.根据二维判据得到的区域液化临界值Y₀=0.236，即当Y大于Y₀，场地为液化，反之不液化.

由公式(1)，峰值加速度PGA与复合地形指数CTI的系数都为正，说明二者越大液化发生可能性越大；平均剪切波速V_s30和场地到河流的距离DR的系数都为负，即二者值越大液化越不容易发生.这与目前认识一致，间接证明了该模型的正确性.

2.2 模型可靠性分析

对模型和公式(1)的可靠性检验，包括整体检验、回归系数显著性检验以及模型系数正确性检验.

本文采用著名的Akaike信息标准(AIC)进行检验，其表达为

(2)

其中K为模型中自变量的数目，S为反应变量类别总数减1，n为观测数量，是所设模型的估计最大似然值的自然对数，其值较大表示拟合较好.- 的值域为0至∞，其值越小说明拟合越好.在其他条件相同时，一个模型的AIC越小说明模型拟合越好.

对于回归系数的显著性检验，本文采用Wald检验，其表达为

(3)

其中S.E.为回归系数估计量的标准差.Wald越大，说明模型回归系数的越显著，系数选取得越好.

模型的整体检验(AIC检验)、回归系数显著性检验(Wald检验)结果如表 1所示.由表 1可知，本文模型各个参数估计的偏差都很小，Wald检验值以及AIC检测值可以表明该模型整体检验较好，回归系数显著.

对于模型的正确性，本文采用ROC(Receiver Operating Characteristic Curve)曲线进行验证，本文模型ROC曲线如图 11和表 2所示.ROC曲线下方的面积(AUC)为0.767，即模型精度为76.7%，表明公式(1)判别液化的区域与实测区域吻合度较高.

2.3 与实测结果对比

将得到的公式代入各个点，可得到各个格网单元评估值，最终结果如图 12，与图 2实际液化调查结果进行对比，可以从直观的角度验证该模型的正确性.具体回判成功率如表 2所示，公式(1)总体回判成功率约为77%，与上述AIC、Wald和ROC的检验结果吻合；液化和非液化回判成功率分别约为73%和78%，表现均衡.

经过分析，模型偏差来源可能有以下几点：

(1) 1976年唐山地震震后三天液化灾害情况解译图可能有偏差.

(2) 选取的建模参数信息精度有待提高，这取决于相关研究成果的更新，分辨率有待提高.

(3) 有些区域先前未液化，但被判别为液化，可能是液化潜在区，只是当时没有迹象，在未来可能会发生液化.这些都需要随着今后各种数据精度的更新进行进一步改进.

3 2011年新西兰地震对本文公式的检验 3.1 2011年新西兰Christchurch地震及其液化情况

2010—2011年间，新西兰坎特伯雷系列地震发生了大范围的液化，该系列地震包括4个重要地震，震级范围5.9~7.1，靠近城市边界及内部的发震断层引发了贯穿Christchurch市的强烈地震动，对建筑物和基础设施造成了极大的破坏，经济财产损失近300亿美元.

该系列地震中，2011年2月22日发生的基督城(Christchurch)M_W6.3地震最具有破坏性，大面积的砂土液化是本次破坏的主因.基督城地区人口38000人，面积约450 km²，近60000栋住宅和建筑遭到破坏，其中20000座建筑破坏严重，桥梁、堤坝、地下生命线设施也受到严重破坏.由于液化破坏严重，难以修复，新西兰政府决定永久放弃部分城区，将近8000座住宅被遗弃.

基督城地震中，市区60%发生液化，液化面积约450 km²，震后调查给出了较为详细的液化区分布，地质资料也可在网上公开获得，为检验本文模型和公式提供了良好条件.

3.2 参数提取

基督城位于Canterbury平原地区，由多条源于南Alps山脉的河流冲积而成的冲积扇形平原.Christchurch市坐落在Waimakiriri古河道上，Avon和Heathcote两条河流穿城而过.Christchurch市及其郊区的大部分浅层土是Waimakariri河沿岸洪水泛滥形成的沉积物，近地表地层具有层状特点，系松散的全新世粉土、砂土等.

本文验证区域面积约500 km²，CTI、V_s30、DR和PGA四个参数的提取结果如图 13所示.

3.3 检验结果

本文提出的区域土壤评估公式(1)对该区域液化预测结果与实测结果的对比如图 14所示.验证区域中，液化场地判别成功率为82.4%，非液化场地判别成功率为84.7%，总体为83.8%，充分证明了本文提出的区域土壤地震液化预测简化方法的可靠性.

4 总结与讨论

发展利用少量简单地质参数给出区域土壤液化的评估方法，震前可预估各地区土壤液化可能性，震后可快速评估震区土壤液化情况，对预防减轻地震灾害以及地震快速响应救援都是有意义的.本文在该方面取得以下进展：

(1) 选取场地平均剪切波速、复合地形指数、场地到河流的距离以及地面峰值加速度这四个参数，分别代表土壤密实程度、饱水状态、地质年代和遭遇的地震作用大小，其中复合地形指数为本文特色，克服了美国NEHRP方法依赖详尽水文地质资料的缺点.

(2) 以我国唐山大地震液化场地为背景，对液化及其附近约9000 km²区域，以震后调查和公开的地质资料和地震动数据为依据，运用GIS获取了近9万个样本点，形成了具备较为典型地质条件的四个参数的数据库.

(3) 运用经典二元Logistic回归方法，建立了区域土壤液化预测模型和评估公式，回归系数AIC检验和Wald检验结果良好，公式对实测结果总体成功率77%，非液化和液化场地成功率分别为73%和78%，表明参数选择正确，模型和公式可靠.

(4) 将建立的公式用于2011年新西兰地震基督城及其郊区近500 km²区域液化评估，场地液化和非液化成功率分别为82%和88%，总体成功率为84%，表明本文方法可行，公式可信.