利用BP神经网络技术建立自贡地区地震动地形效应模型及其在汶川地震中的效验

周红， 常莹

1 中国地震局地球物理研究所， 北京 100081 2 矿冶科技集团有限公司矿山工程研究设计所， 北京 100160

0 引言

20世纪地震中地形效应所产生的严重危害，就已被观测和记录到：如1971年的San Fernando地震(Trifunac and Hudson，1971)，1985年的Chile地震(Çelebi, 1987)，2003年的土耳其Bingol地震(Akkar et al.，2005)等，震害调查结果均显示位于山顶和山脊上的建筑，比位于山体底部的建筑遭受到更严重的损害.研究地形效应对于地震工程学及地震学发展非常重要.我国抗震设计规范(GB50011-2010)简单的把盆地和地形效应最大放大系数定为1.6，欧洲设计规范 (EuroCode 8, CEN 2003)及意大利设计规范(NTC-08 2008)基于理论计算结果规定场地最大放大系数为1.4，对于不同的地区不同的城市都采用一个统一的放大系数，看似简单，实则与实际情况相悖.

由于观测记录的缺乏，从台站记录资料中直接研究地形效应，得到一个科学的普遍性结论非常困难(Wang et al.，2004).因此利用数值模拟技术，开展地形效应形成机制和时空变化规律的研究，目前仍是该研究课题的主要方法.此项研究开始于20世纪70年代，由于计算理论和计算机能力的限制，其研究主要围绕二维地形效应的效应问题进行(Aki and Larner，1970；Sánchez-Sesma and Rosenblueth，1979； Wong and Jennings，1975； Wong，1982；Dravinski，1982；Campillo and Bouchon，1985).Boore等(1981)、Bouckovalas和Papadimitriou(2005)研究了台阶地形坡面上下拐角对地震动的散射作用；Ashford和Sitar(1997)研究了二维坡顶对地震动的响应；Bard(1982)研究了高斯型山脊对SH波的响应；Bouchon等(1996)、廖振鹏等(1981)研究了孤立山地的放大作用；周红和陈晓非(2007)研究了二维山和谷地形对于Rayleigh波的影响；Geli等(1988)总结了前人对二维孤立山体的研究，讨论了不同山体的高度h和山体的基底的宽度w之比对地震动的影响，认为地形尺度和地震波波长是地形效应的重要影响因素，相似的结论在其他文献里也被证实(Ashford et al.，1997；周红和陈晓非，2007).近几年实际区域内三维地形效应的研究得到了开展：Ma等(2007)研究了San Gabriel山脉长周期的地形效应；Lee等(2008)研究了台北周边长周期(约0.2 Hz)的地形效应；Lee 等(2009)结合高分辨率地形资料，研究了台湾阳明山地区不同地震参数对高频(10 Hz)地震动地形效应的影响；蒋涵等(2015a,b)研究了芦山地区山脊和坡度与地形效应的关系.

现在已有的场地效应预测模型都是单参数、体现平均效果的统计模型，如林柏伸和李锡踶(2003)建立了地形比高因子与地形效应之间的平均放大系数的统计模型；蒋涵等(2015a；2015b)建立了芦山地区地形坡度与地形效应之间的统计关系的数学模型；Kim等(2018)建立了三个独立关系模型：地形效应分别与地形坡度、地形曲率、地形高程的统计关系；Maufroy等(2015)考察了不同频率的场地效应与地形曲率的统计关系.地形构造是非常复杂的，不单单与地形高程、地形坡度等简单参量有关；复杂的地形效应亦不仅仅与简单的单一参量有关，而是由复杂的参数体系共同决定，因此依据单参量获得的地形效应统计预测模型难以准确代表实际地形场地效应.本论文将不再使用简单的统计方法获得平均的预测模型，而是引入神经网络技术，得到地形效应每个场地的一对一的预测模型.

国家“十五”重大项目《中国数字地震观测网络项目》在四川盆地南部的自贡市自流井区西山公园内布置了地形影响台阵，目的是为了研究地形对地震动的影响，从2007年开始运行.2008年汶川地震中，该台阵记录了良好的主震加速度记录，提供了宝贵的研究地形效应的实际资料.基于这些观测记录学者们已获得了一些研究成果，如基于分析地震动记录参数特征，研究地震动放大效应(阿布都瓦里斯等，2013；王伟等，2015，2016；王海云和谢礼立，2010)；基于二维数值模拟，研究地震动放大效应(杨宇等，2011；张建毅等，2012；唐晖等，2012).

本文针对自贡台阵区域，利用三维谱元数值模拟技术，获得研究区内的地表地震动，计算相应的地形放大系数，接着引入BP神经网络技术，建立地形效应预测模型，并将地形效应预测结果与自贡台阵实际记录的加速度峰值(PGA)放大结果相对比，给出了差异范围和可能的原因.

1 方法介绍

本文首先采用谱元数值模拟技术(Komatitsch and Vilotte，1998；Komatitsch et al., 2000；Komatitsch and Tromp，2002；Zhou and Chen，2010；Zhou and Jiang，2015)计算研究区域内的地面运动场，进而计算得到地表地震动地形效应放大系数；并引入BP神经网络非线性建模技术，建立研究区内各点地形特征与地震动放大比值的关系，完成地形效应多参数模型的建立.下面分别简介研究中涉及的谱元模拟方法和BP神经网络建模方法.

1.1 地震波谱元模拟方法

在地下介质中传播的地震波，其位移场u满足如下的波动方程：

(1)

谱元法是一种微分方程变积分方程的数值方法，对方程(1)点乘任意试错函数w，并采用分步积分可以得到方程(2)：

(2)

这里Ω是物理的计算域，Γ是它的边界，包括自由地表边界Γf和人工边界Γabs，由于地表应力为零，Γf在分部过程消失，因此上式自然满足地表应力条件.

对于三维问题，为了完成方程(2)的积分，需要将整个计算区域Ω剖分离散成非重叠的不规则形状的六面体小单元，在计算过程中需要将不规则形状的剖分单元映射成规则单元，每个单元采用Gauss-Lobatto-Legendre(GLL)点进行积分计算，结合GLL积分规则，可以得到关于计算地震波位移的矩阵方程组：

(3)

其中M是质量矩阵，K是刚度矩阵，F是源项，U表示未知量位移的向量.GLL积分规则的引入，使得M变为对角线矩阵，这种对角性使得谱元方法(SEM)易于开展并行计算.

1.2 BP神经网络建模方法简介

BP(Back Propagation)神经网络是一种最常用的人工神经网络技术，它是Rumelhart和McCelland(Rumelhart et al., 1986; Anderson and Rosenfeld, 1988)为首的科学家小组提出来的一种按误差逆传播算法的前馈网络，该网络可以学习和存贮大量的输入、输出映射关系，而无需预先已知映射的数学方程.其学习规则是使用最速下降法等搜索技术，通过反向传播不断调整网络的权值和阈值，使网络的误差平方和达到最小.BP神经网络模型通常具有三层结构：输入层(input)、隐层(hide layer)和输出层(output layer)，其计算流程包括信号的正向传播及误差的反向传播.本论文利用BP神经网络，建立自贡西山公园平面波入射条件下的多参数地形效应.

2 自贡西山公园地震动地形效应的BP神经网络模型的建立

2.1 自贡西山公园地形地貌

国家“十五”重大项目《中国数字地震观测网络项目》在四川盆地南部的自贡市自流井区西山公园内布置了地形影响台阵，目的是为了研究地表地形对地震动的影响.台阵在水平自由地表场地、山脚、山底及山脊等位置布设了8个观测点，其中1号测点为水平自由地表土层测点，2号测点为山脚基岩测点，其余测点为山体基岩测点，各测点位置如图1所示.

图1 研究区域地形及台阵分布圆圈及标号标示了1～8号观测点位置.测线A和B分别沿着纬度线29.335°N和29.344°N，用于标示图3展示的两组100 m间距测线点所沿测线位置.Fig.1 The study area and the locations of seismic stations in the Zigong arrayBlue circles and numerical labels show the positions of the eight stations. The profiles A and B, along the latitudes of 29.335°N and 29.344°N, respectively, show the 100-m-interval sites of the simulated time series in Fig.3.

为了利用1～8号台站记录研究地形对地震动的影响，我们选择包含西山公园在内的2.16 km×2.88 km区域作为研究区域，建立3D地形计算模型(如图2a)，很显然，研究区内地形起伏较小，最高海拔与最低海拔之差小于100 m，因此研究区属于地势较为平坦的区域.

为了去除地表非均匀介质对地形效应的干扰，我们基于玉皇观 (余嘉顺等, 2017)的勘探结果设定3D计算模型的物性为S波速度600 m·s-1单一的基岩介质；同时忽略不同方向传播特性对地震动的影响，设定垂直入射的S波为激发源，其时间函数取为简单的Ricker函数，主频为10 Hz，频宽0～30 Hz(图2b).

图2 计算模型(a) 计算介质模型; (b) 归一化的射入波频谱.Fig.2 The computational model(a) The calculated medium model; (b) The normalized incident wave spectrum.

2.2 地震动模拟结果的地形效应

图3是图1中剖面A、B对应的Ricker子波激发下的地表地震动模拟时程.在地形高程中(图3A，图3B)的地形， 显示地形存在多处峰谷，但是高程差异较小，最大不超过60 m，对应的地表地震动(图3A1)的到时、波形差异都很小，可以说没有明显差异；尽管局部峰谷起伏很小，却存在一个较长尺度的隆起，峰谷高程差接近100 m，对应的地表地震动(图3B1)到时差异明显，各台站波形也不尽相同.我们将利用这些模拟的地表地震动来研究地形对地震动的放大效应.

图3 测线A、B(所在位置见图1)地表地形与模拟的地表加速度时程(A1) 沿测线A的地表加速度时程； (A2) 测线A高程； (B1) 沿测线B的地表加速度时程; (B2) 测线B高程.Fig.3 The time series of the simulated surface ground motions and the topographic elevation along the profiles A and B in Fig.1.(A1) is the simulated ground accelerations along the profile A; (A2) is the topographic height along the profile A; (B1) is the simulated ground accelerations along the profile B; (B2) is the topographic height along the profile B.

图4为计算模型(图2a)在垂直入射Ricker横波作用下的地形对地表地震动峰值加速度PGA的放大系数空间分布.地形对地震动放大系数的计算目前没有统一的计算标准，本文中的放大系数采用方程(4)计算，

(4)

其中V为研究区内各计算场点的PGA或者频谱值， ref(V)为参考点的PGA或者频谱值.我们选择图1中位于山脚的2号基岩测点位置为参考点，参考量的选取不改变放大比与地形的关系，但影响地形放大比的绝对值.

图4 地形对地震动的PGA放大系数的空间分布Fig.4 The spatial distribution of topographic magnification factors of PGA

以PGA为例进行讨论，图4为PGA放大比分布，图中可见较大的放大系数(红色表示)分布比较零星，多位于地形起伏的局部隆起位置上，小于1的地形放大系数分布范围较广，多位于坡面或者坡底.图5是图1中的两条剖面A、B对应的地形及其模拟地震动计算得到的PGA放大系数，可以看到地形海拔高度与地形放大系数成复杂的相关性，每个局部地形的高点几乎都对应着放大系数峰值，地形的局部低点对应地形放大系数低值谷点，但是地形的区域最高峰点，并不对应最大放大系数，同样地形的区域最低点，也不对应放大系数的最小点，地形与放大系数具有相关性是确定的，但是相关性又难以直接表述.

图5 剖面A、B对应的地形及其模拟地震动PGA地形放大系数Fig.5 The topography along the profiles A and B and the topographic PGA magnification factors of simulated ground motions

图6a显示了计算区域地形高程与PGA放大系数的对比，对于某一个高程，放大系数从小于1到大于1都存在，而这种现象存在于本研究区域的所有不同的高程上，并没有哪个高程占优.对于地形坡度与放大系数的关系上(图6b)，也存在这种现象，图6显示似乎地形放大系数与地形高程、地形坡度没有确定的关系，或是仅是随机的关系，难以通过地形高程、坡度预测地形的放大系数，事实真是如此吗？

图6 计算区域地形高程、地形坡度与放大系数的对比(a) 地形高程与放大系数的对比； (b) 地形坡度与放大系数的对比.Fig.6 Comparisons of the regional topographic elevation and slope with the topographic magnification factor(a) Comparison between the topographic elevation and magnification factor; (b) Comparison between the topographic slope and magnification factor.

2.3 地震动地形效应的BP神经网络预测模型的建立

图4—图5显示PGA的地形效应与山脊、次级山脊、谷、次级谷等地形构造有明显的正相关性，但是图6显示地形效应与地形高程、地形坡度没有明确的关系，任何高程，任何坡度的地形都有可能对地震动放大，也有可能对地震动缩小.本节将利用BP神经网络技术，通过猜测变量的方式，找到何种地形特征是影响地形对地震动放大效应的主要因素，并且建立研究区域地形效应的定量预测模型.

建立BP神经网络模型，需要首先确定神经网络模型的参量，前人的研究成果认为地形效应与地形高程、坡度、地形尺度、地形宽高比的变化等相关，对于真实的地形，除了地形高程、坡度，其他参量如地形尺度、地形宽高比都无法用特定的量描述，那么用什么参量可以比较充分的表达地形的各种特征呢？Zhou等(2020)的研究成果显示：以地形高程、地形南北向和东西向的梯度作为基本参量，以地形南北向和东西向的高阶梯度作为地形尺度相关的参量，建立地形效应模型有效.这里我们继续采用这些参量.其中地形南北向和东西向的多阶梯度采用数值微分法计算(刘敏等，2007).该方法利用计算点周围的8个点所组成的窗口，即3×3的窗口来计算地形多阶坡度和梯度(图7).

图7 坡度计算示意图Fig.7 The schematic diagram of slope calculation

每个窗口中心为被计算点，该点的坡度的定义式如下：

(5)

其中fx和fy为南北向和东西向的梯度，可以通过方程(6)获得，即：

(6)

hi为计算窗内各点高程，g为单元的长度.南北向梯度fx和东西向梯度fy，地形的二阶梯度也采用同样的手段计算，只是把高程h替换为fx和fy.

选择三层BP模型，设定学习速率为0.05，最大训练轮回为2000；学习均方误差为10-4搭建神经网络进行模型训练.地形不仅放大地震动的PGA，也影响着地震动的各个频率成分，地震动加速度反应谱(PSA)能够反映地震动不同频率成分产生灾害的强弱，因此被广泛采用.本文选择三组常用频率的PSA(0.3 Hz、1.0 Hz、3.0 Hz)和PGA进行建模实验，其中PGA相当于无限频率对应的反应谱.建模过程首先分别将PGA和PSA放大系数随机地分成两组样品：随机抽取80%的PGA和PSA放大系数定为训练样品，而剩余的20%的PGA和PSA放大系数被定为预测样品.图8为不同参量建立的神经网络BP模型预测结果(红线)与PGA、PSA放大系数(黑线)的对比，这里纵轴表示PGA、PSA放大系数，横轴是随机取样的地表场地点.图8(a—i)(i=1，2，3，4)是以地形高程、地形坡度(h,slope)作为参量BP模型预测结果与由方程(4)计算的放大系数的对比.以(h,slope)为参数BP模型预测的放大系数近似等于1，与实际地形放大系数相关性差，表明地形效应与地形高程、地形坡度没有明确的关系，或者仅是随机的关系，与图6的结论一致.图8(b—i)是以地形高程、地形南北梯度、地形东西梯度(h,fx,fy)作为参量的预测结果与实际地形效应对比，预测效果较图8(a—i)有所提高，但是预测结果依然近似等于1，这点再次映证了图6的结论.图8(c—i)是以地形高程、地形南北梯度、地形东西梯度、地形南北二阶梯度、地形东西二阶梯度(h,fx,fy,fxx,fxy,fyx,fyy)作为参量的BP预测结果与实际地形放大系数对比，与图8(a—i)和8(b—i)相比较，BP预测结果有了大的改善，预测的地形效应值已非常接近由方程(4)计算的实际地形放大值，这也表明地形南北二阶梯度、东西二阶梯度所反映的小尺度地形特征在地形放大效应中起了较大的作用.图8(d—i)是以地形高程、地形南北梯度、地形东西梯度、地形南北二阶梯度、地形东西二阶梯度、地形南北三阶梯度、地形东西三阶梯度(h，fx,fy,fxx,fxy,fyx,fyy,fxxx,fxxy,fxyx,fxyy，fyxx，fyxy，fyyx，fyyy)为参量的BP模型预测结果与实际放大系数对比，该组参量的预测结果和图8(c—i)都接近直接计算的放大值，说明加入更小尺度的地形特征，可以提高预测结果的准确性.图8(a—i)到图8(d—i)预测值的变化形态显示：对地形效应影响较大的成分隐含在地形南北二阶梯度、东西二阶梯度中，而其他尺度地形特征的梯度，如一阶梯度、三阶梯度，影响作用较小.均方根RMS的差异的定量比较(表1)更清楚地说明这点，图8(b—i)到图8(c—i)RMS的减少量是所有模型里最大的，而图8(c—i)—图8(d—i)RMS的减少量比图8(b—i)—图8(c—i)RMS的减少量小很多，因此对于该计算地区，地形更小尺度的三阶地形梯度对于预测模型只有较小的、有限的影响.考虑到计算成本的因素，作者认为，针对本研究区域，取地形高程、地形南北梯度、东西梯度、南北二阶梯度、东西二阶梯度(h,fx,fy,fxx,fxy,fyx,fyy)作为参量建立BP神经网络模型，可以达到该研究区地形放大系数建模的需求.

地震动的放大效应不但与地形相关，与地震波传播方向也密切有关，入射方向的改变会影响地形效应放大比值的大小，但不会改变地形南北梯度、东西梯度、南北二阶梯度、东西二阶梯度对地形效应的作用模式，附图A显示了这点.附图A为45°入射S波各频率不同参数系BP模型预测的结果，与图8各点一致.相较于图8放大比值存在差异，但地形效应与地形参量的定量关系模式继续有效.关于更多不同角度入射激发的地形效应的特殊规律，我们将在专门的论文里做详细讨论，这里不做展开.

表1 不同参数预测模型的RMS对比Table 1 RMS comparison of the difference prediction models

图8 不同参量神经网络模型结果与实际计算的四组地形放大系数(0.3 Hz,1.0 Hz,3.0 Hz PSA和PGA) 的对比Fig.8 Comparisons of magnification factors predicted from different parametric neural network models with the topographic magnification factors of PSA at 0.3 Hz, 1.0 Hz PSA, 3.0 Hz and PGA calculated from simulated ground motions

2.4 BP神经网络模型预测结果与自贡台阵观测记录地形效应的对比

图9为自贡地形台阵2～8号台站的东西向和南北方向加速度观测记录，台站2的记录持时120 s，PGA约为25 cm·s-2，频谱宽度0～45 Hz.为了与模拟地震动对比，观测记录采用Butterworth 30Hz低通滤波.

图9 自贡2～8号台EW向和NS方向加速度记录Fig.9 Acceleration records in the EW and NS directions of Zigong station No.2～8

分别以谱元法数值计算和汶川地震观测记录的2号台站(st2)PGA为方程(4)中放大系数的参考点数据，通过图8c BP神经网络模型预测的台站st3-st8地形放大系数，与汶川地震观测记录获得的st3-st8地形效应放大系数做比较.台站st1为土层，在不考虑场地效应下不参与地形效应对比.表2给出了对比结果：预测模型和观测记录均以st2场点作为参考点，对于NS方向地形放大值相对差最大的为台站st8，相对差47%；EW方向地形放大系数相对差最大的在台站st7，相对差48%；相对差最小是台站st5，NS方向相对差0.8%，EW方向没有差异，其他台站EW和NS方向预测与实际相对差异在10%～48%之间.从绝对差异上看，预测的放大系数普遍比实际值小.

表2 台站3～8的观测PGA放大系数与BP预测的对比Table 2 Comparison between the actual magnification factors of PGA of stations 3～8 and the theoretical BP predicted coefficients

引起这些差异的原因可能为如下四个：

1)地形高程模型分辨率不足，我们建立的地形计算模型，资料来源于STRM1(30 m分辨率的Shuttle Radar Topography Mission资料)，最小分辨率30 m，影响实际记录的地形放大值可能是小于分辨率尺度的地形特征，我们采用的高程资料不能体现，其对地震动的影响也没在模拟结果中出现；

2)模拟记录由简单的Ricker子波产生，本文的模拟记录是垂直入射的Ricker子波的结果，模拟的地震动比实际观测地震动简单很多，尽管将对比的实际记录滤波到0～30 Hz，只能保证对比的频宽近似一致，不能保证时程上变化一致，也不能保证PGA相对一致，因此基于模拟结果建立的地形效应BP预测模型，显然有别于实际地形效应；

3)地表介质非均匀性，实际记录是在非均匀性的地表介质下产生的，本研究区域为单一介质；

4)传播方向的影响，BP模型是对垂直入射地震动地形效应的预测，实际记录地形效应可能包含传播方向的影响.

尽管如此，本文建立的预测模型给出的地形放大的结果，给出了对实际地震地形效应预测误差的范围，相信随着更多场地预测模型的建立，更多复杂地下构造模型、震源模拟的使用，会得到更为贴合实际的预测模型，为真正的地形效应预测提供更为可靠的结果.

3 结论

本文旨在确认自贡地区地形效应的决定因素，建立其定量预测模型.我们首先建立研究区域单一地形模型，模拟垂直入射平面S波的地表地震动，获得地形对地震动放大系数的空间分布，空间分布显示地形放大与山脊、次级山脊关系密切，但地形高程、地形坡度与地形效应又呈现出无关的矛盾现象.为此我们引入BP神经网络非线性建模方法，通过对比分析确定了地形效应的主要影响因素，即地形效应与地形高程、地形的南北向、东西向梯度及其高阶梯度密切相关，其中地形高程代表大尺度信息，地形的南北向、东西向梯度表示小尺度地形信息，而地形的南北向、东西向高阶梯度包含了更小尺度的地形信息.预测结果的数值对比还显示，对本研究的模型而言，这些参量中起最主要作用的是南北向、东西向二阶梯度，它的加入明显的区分开各个场地的地震动放大系数.论文最后利用地形高程、地形的南北向、东西向梯度及其高阶梯度参量建立的研究区域的BP神经网络建模预测了自贡地形台阵的地形放大系数，并与汶川地震获得的实际NS、EW向地形放大系数做了对比(附图A)，总体来说预测结果相对较小，台站st5地形预测结果与实际最为符合，差异约为1%；而其他台站预测结果和实际地形放大系数相对差异在10%到48%之间，这种差异可能来自1)地形数据的分辨率不够高；2)入射波形简单；3)地下介质模型与实际介质差异大；4)忽略了地震波传播方向的影响.

总结来看：地形效应不仅与地形坡度有关，它受多种地形变化特征的影响，因此地形效应预测模型需要由多参数来表达，而非单一参数模型可以描述.对于自贡地区的Str1分辨率地形，这组参数由地形高程、地形的南北向、东西向梯度及其二阶梯度构成，以此建立的研究区地形效应BP神经网络模型，预测的地震动放大系数与实际的地形效应相对差异在0～48%之间.未来，随着高分辨率地形资料使用，并且考虑多方向、复杂入射波以及非均匀地下介质，地震动地形效应预测模型将得以改善，关于这些因素的影响将在以后的研究中逐一讨论.

致谢中国地震局工程力学研究所“国家强震动台网中心”为本研究提供数据支持；国家超级计算天津中心和广东中心提供超算计算平台.

附录A

地震动的放大效应有很多影响因素，比如地形、地球内部的速度横向非均匀性，以及地震波入射方向等，地形起伏本身对地形效应的影响，对于不同的入射波、不同的入射波影响模式保持一致，地形南北梯度、东西梯度、南北二阶梯度、东西二阶梯度是地形效应建模的重要参量.

附图A 不同地形参量神经网络模型预测结果与实际计算的四组地形放大系数的对比order=0表示地形高程h为建模参量；order=1表示(h，fx,fy) 为建模参量； order=2表示(h,fx,fy,fxx,fxy,fyx,fyy)为建模参量；order=3表示(h,fx,fy,fxx,fxy,fyx,fyy,fxxx,fxxy,fxyx,fxyy,fyxx,fyxy,fyyx,fyyy)为建模参量.Appendix Fig.A Comparisons of magnification factors from different parametric neural network models with the actual topographic magnification factors.order=0,1,2,3 respectively denote that the established BP models depend four sets of parameters h,(h,fx,fy),(h,fx,fy,fxx,fxy,fyx,fyy) and (h,fx,fy,fxx,fxy,fyx,fyy,fxxx,fxxy,fxyx,fxyy,fyxx,fyxy,fyyx,fyyy).