Study on 2D in-plane HVSR simulation and application with transverse inhomogeneous body scattering
-
摘要: 为分析场地的横向不均匀性对水平竖向谱比HVSR曲线产生的显著影响，本文基于Sánchez-Sesma等提出的扩散场方法，通过计算总波场格林函数虚部对二维沉积地形上的HVSR曲线进行模拟。格林函数虚部则通过刚度矩阵和平面内斜线格林函数采用间接边界元方法进行求解。对二维沉积地形和相应的一维层状半空间的HVSR曲线进行了参数分析，着重讨论了沉积地形的形状、相对计算点位置等因素对HVSR曲线的影响规律。结果表明：沉积地形内外材料阻抗比对HVSR曲线的影响最为显著；随着沉积地形内外材料阻抗差异和沉积侧界面坡度的增大，沉积地形上HVSR曲线的第一峰值点的频率显著增大至相应层状半空间结果的3.3倍，同时HVSR曲线的形态呈现出平台现象；随着计算点到沉积边界距离的减小，HVSR曲线高频段幅值相对较大。根据本文得到的局部地形对HVSR曲线的影响规律，在进行场地勘探时可采用HVSR方法初步确定局部地形的分布位置以降低勘探成本。Abstract: In order to analyze the significant influence of lateral inhomogeneity of site on horizontal-to-vertical spectral ratio (HVSR) curves, the diffuse field approach proposed by Sánchez-Sesma et alwas adopted to simulate the HVSR curves of 2-D sediment topography by calculating the imaginary part of Green’s functions of total wave field. The imaginary part of Green’s functions was solved by the dynamic stiffness matrix and in-plane inclined Green’s functions based on the indirect boundary element method (IBEM). The HVSR curves of 2-D sediment topographies and corresponding 1-D layered half-space were compared, the influences of sediment topography shapes and the relative position of calculation points on the HVSR curve were discussed in detail. The results show that the effect of impedance ratio between inside and outside materials of sediment topography on HVSR is the most significant; With the increase of the impedance ratio and the slopes of the interface on the sediment side, the frequencies of the first peak of HVSR curves increase significantly, which can be up to 3.3 times of the corresponding layered half-space results, simultaneously, platform emerges on HVSR curves; Amplitudes of HVSR curves in high frequency band increase with the decrease of distances from the calculation points to the sediment boundary. According to the results obtained in this study, the HVSR method can be used to preliminarily determine the place where local sediment topography exists. From this aspect, the cost of regional geophysical investigation can be reduced visibly via HVSR method.
-
引言
场地背景噪声和地震动与地球表面地层结构特性有着紧密的关联。地表运动混合了各种形态的入射波和散射波，同时包含了波源、传播路径和地下结构等信息。对背景噪声或地震动记录采用水平竖向谱比（horizontal-to-vertical spectral ratio，缩写为HVSR）方法进行场地效应研究，是当前国际上的研究热点。HVSR曲线在场地特征参数测定、地震区划、场地响应评估、场地速度结构反演及地下采空区和岩溶发育地区探测等工作中具有重要的理论意义和应用价值（温瑞智等，2015；荣棉水等，2016，2018；Piña-Flores et al，2017）。
目前HVSR方法的应用研究大多基于一维层状半空间场地的假设，但在实践中，局部地形的存在会对HVSR曲线的形态和峰值频率造成显著影响，干扰HVSR曲线的应用。王伟君等（2011）指出，受到盆地几何形状（尤其是比较狭窄盆地）的二维或三维效应的影响， HVSR曲线峰值频率会产生偏移，进而影响采用HVSR方法得到的土层卓越频率。Guéguen等（2007）研究表明，对于狭窄的法国格勒诺布尔沉积盆地，较小的盆地宽度和沉积厚度之比使得二维和三维效应明显，某些情况下用峰值频率得到的沉积层厚度误差可达10%，而在盆地边缘误差最高可至50%。此外，Stanko等（2017）对位于克罗地亚北部沉积盆地上的瓦拉日丁城进行了微震观测，指出断层界面对HVSR曲线有方向极化作用。上述研究表明，对于横向不均匀情况明显的场地，不宜再使用一维层状场地模型进行HVSR曲线的计算和分析。
同时，HVSR方法在探测局部地形和地下结构方面也表现出了巨大的潜力， Chavez-Gárcía等（1997）指出可以采用HVSR方法评估局部地形场地效应。即根据测点采用HVSR方法确定的局部地形场地效应的强弱，判断测点与局部地形的相对位置和局部地形特征的相关信息。若能对局部地形中的HVSR曲线进行准确模拟，并在此基础上细化分析出局部地形因素的影响规律，将HVSR方法进一步应用于局部地形场地的探测和反演，对实际工作具有重要的理论意义和价值。
目前，针对局部地形的HVSR研究大多通过场地的实测地震记录进行计算分析。研究人员基于实测数据分析了局部地形对HVSR曲线的影响规律，并得出了相应结论（Gosar，2007；卢育霞等，2017；Napolitano et al，2018；林国良等，2019）。但基于实测数据得到的HVSR曲线只针对测点所在位置和地形，难以得到一般性的结论。为了总结局部地形对HVSR曲线的影响规律，需要较准确的数值模拟结果作为辅助参考，细化分析横向不均匀体的各个因素对HVSR曲线影响的敏感性和影响规律。部分研究人员利用数值方法对局部地形上的HVSR曲线进行了模拟（Chavez-Gárcía et al，1997；Uebayashi et al，2012）。但由于场地中背景噪声源和地震震源的复杂性，数值方法难以准确模拟场地所受到的激励，得到的HVSR曲线往往与实测结果相差较大，很难真实地反映局部地形对HVSR曲线的影响。
Kawase等（2011）和Sánchez-Sesma等（2011）通过引入扩散场理论解决了入射波源的复杂性问题，在场地波场具有扩散场性质的基础上，间接得到了HVSR曲线，提出了模拟一维层状半空间上HVSR曲线的扩散场法。此后，Matsushima等（2014）根据扩散场理论采用谱元法计算了三维沉积地形上的HVSR曲线并分析了局部地形使HVSR曲线产生的方向性特征。Perton等（2018）采用扩散场法对具有横向不均匀性火山口地形的HVSR曲线进行了模拟，并对地层速度进行了反演。但以上均未进行局部地形对HVSR曲线的影响规律研究分析。
为了更为准确、快速地模拟局部地形场地上的HVSR曲线，本文在扩散场理论的基础上运用计算精度高、速度快的间接边界元方法（indirect boundary element method，缩写为IBEM）对沉积地形上HVSR曲线进行模拟，总结出沉积地形特征对HVSR曲线的影响规律。在此基础上，探讨了根据HVSR曲线特征反向确定局部地形信息的可能性，为采用HVSR方法勘探局部地形和利用HVSR曲线反演二维地下结构等潜在应用提供基础。
1. 扩散场法计算原理
日本学者Nogoshi和Igarashi （1971）提出了根据地震作用下的场地响应确定场地特征频率的方法。Nakamura （1989）在此基础上总结出了HVSR方法并将其进行了推广。在推导HVSR方法的过程中，Nakamura提出了两个假定：① 覆盖土层对基岩处的竖向位移分量不产生放大效应；② 基岩处的水平和竖向位移分量相等。因此，HVSR曲线的幅值等价于地表水平位移与基岩处水平位移的比值。相应的，HVSR曲线的峰值频率等于单一层状半空间场地一维速度结构的卓越频率，即
$$ {f_n} {\text{＝}} n{}{\frac{v_{\rm{S}}}{4h}} {}{\text{，}} $$ (1) 式中，vS和h分别为土层的剪切波速和厚度，n为覆盖土层的第n阶模态，取为1，3，5，···。然而Nakamura （1989）提出的假设在场地中存在局部地形时往往不能成立，这是由于局部地形产生的散射波增加了场地中波场的复杂性，实际情况会与假设不符（王伟君等，2011）。
为了将HVSR方法理论扩展到存在局部地形的一般场地，Kawase等（2011）和Sánchez-Sesma等（2011）提出了扩散场法并为HVSR曲线给出了新的解释。在扩散场法中，场地的背景噪声波场和地震动波场被认为是具有能量均分性质的扩散场，水平与竖向位移比值代表水平方向与竖直方向能量密度的比值。在一个可被视为扩散场的局部地形场地的波场中，可采用扩散场法对HVSR曲线进行计算。
扩散场的概念最早是在声学领域中提出的，考虑弹性介质区域受到一组不相关的随机力，所得到的波场相关性质与多次散射波场是等效的（Sánchez-Sesma，2017）。由于后者可以用类扩散方程很好地描述，Sánchez-Sesma等（2011）用扩散场这一概念描述地层中的噪声场。扩散场的表征之一是能量均分，这一概念最先由Weaver （1982）在全空间中提出，即在弹性介质中的各个方向能量相等，各个波形所具有的能量之比为定值。此后Weaver （1985）和Perton等（2009）将其拓展到半空间的情况。在能量均分的前提下，Sánchez-Sesma等通过解析解计算证明格林函数张量虚部与方向能量密度成正比，并给出了用位移互相关表示格林函数张量虚部的表达式（Sánchez-Sesma，Campillo， 2006；Sánchez-Sesma et al，2008），即
$$ {{\rm{Im}}} \left[ {{{{G}}_{{{i}},{{j}}}}{\text{（}}{x_B},{x_{A}},\omega {\text{）}}} \right] {\text{＝}}\frac{{ - \rho {\omega ^2}}}{{4\mu {E_{\rm{S}}}}}\left\langle {{u_{i}}{\text{（}}{x_{A}},\omega ,\theta {\text{）}}u_{j}^*{\text{（}}{x_B},\omega ,\theta {\text{）}}} \right\rangle {\text{，}} $$ (2) 式中：格林函数
${{{G}}_{{i},{j}}}{\text{（}}{x_B}{\text{，}}{x_A}{\text{，}}\omega {\text{）}}$ 表示在点xA沿j方向处施加一频率为ω的单位简谐荷载，在点xB接收到的沿i方向的位移响应；μ和ρ分别表示材料的剪切模量和密度；${u_{i}}{\text{（}}{x_{A}}{\text{，}}\omega {\text{，}}\theta {\text{）}}$ 和${u_{j}}{\text{（}}{x_B}{\text{，}}\omega {\text{，}}\theta {\text{）}}$ 分别表示点xA和点xB受到入射角为θ单一平面波产生的位移响应；ES表示半空间中剪切波具有的能量密度；尖括号表示对所有方向入射波产生的位移互相关取方位平均；上标“*”表示进行复共轭运算。在三维情况下，令式（2）中的接收点xB和作用点xA重合，此时，格林函数虚部代表着单位简谐荷载对介质输入的能量，即$$ {E_m}{\text{（}}{x_A}{\text{）}} {\text{＝}} \rho {\omega ^2}\left\langle {{u_m}{\text{（}}{x_A}{\text{）}}u_m^*{\text{（}}{x_A}{\text{）}}} \right\rangle {\text{＝}} \frac{ - 2\pi \mu {E_{\rm{S}}}}{k} {\text{×}} {{\rm{Im}}} \left[ {{G_{mm}}{\text{（}}{x_A},{x_A}{\text{）}}} \right]{\text{，}} $$ (3) 式中：Em表示方向能量密度（directional energy density，缩写为DED），m＝1，2，3，分别代表平面内水平方向、平面外水平方向和竖直方向；k表示剪切波波数k＝ω/vS。假设微震的波场属于散射场，那么地表点x处的HVSR可以用方向能量密度Em表示为（Arai ，Tokimatsu，2004）
$$ \frac{{\rm{H}}}{{\rm{V}}} {\text{（}}\omega {\text{）}} {\text{＝}} \sqrt {\frac{{{E_1}{\text{（}}x{\text{，}} \omega {\text{）}} {\text{＋}} {E_2}{\text{（}}x{\text{，}} \omega {\text{）}}}}{{{E_3}{\text{（}}x{\text{，}} \omega {\text{）}}}}} \;\;{\text{，}} $$ (4) 将式（3）代入式（4），可以得到用格林函数虚部表示的HVSR，即
$$ \frac{{\rm{H}}}{{\rm{V}}}{\text{（}}\omega {\text{）}} {\text{＝}} \sqrt {\frac{{{{\rm{Im}}} \left[ {{G_{11}}{\text{（}}x{\text{，}} x{\text{，}} \omega {\text{）}}} \right] {\text{＋}} {{\rm{Im}}} \left[ {{G_{22}}{\text{（}}x{\text{，}} x{\text{，}} \omega {\text{）}}} \right]}}{{{{\rm{Im}}} \left[ {{G_{33}}{\text{（}}x{\text{，}} x{\text{，}} \omega {\text{）}}} \right]}}} \;\;{\text{，}} $$ (5) 对于二维平面内情况，不存在平面外水平方向位移u2和能量密度E2，则荷载作用点的能量密度可以写成（Sánchez-Sesma et al，2011）
$$ {E_m}{\text{（}}{x_{A}}{\text{）}} {\text{＝}} \rho {\omega ^2}\left\langle {{u_m}{\text{（}}{x_{A}}{\text{）}}u_m^*{\text{（}}{x_{A}}{\text{）}}} \right\rangle {\text{＝}} - 4\mu {E_{\rm{S}}} {{\rm{Im}}} \left[ {{G_{mm}}{\text{（}}{x_{A}}{\text{，}}{x_{A}}{\text{）}}} \right] {\text{，}} $$ (6) 其中，m＝1和3分别代表平面内水平方向和竖直方向。地表HVSR则可表示为（Matsushima et al，2014）：
$$ \frac{{\rm{H}}}{{\rm{V}}}{\text{（}}\omega {\text{）}} {\text{＝}} \sqrt {\frac{{{E_1}{\text{（}}x{\text{，}}\omega {\text{）}}}}{{{E_3}{\text{（}}x{\text{，}}\omega {\text{）}}}}} {\text{＝}} \sqrt {\frac{{{{\rm{Im}}} \left[ {{G_{11}}{\text{（}}x{\text{，}}x{\text{，}}\omega {\text{）}}} \right]}}{{{{\rm{Im}}} \left[ {{G_{33}}{\text{（}}x{\text{，}}x{\text{，}}\omega {\text{）}}} \right]}}} \;\;{\text{.}}$$ (7) 由式（4）和式（7）可以看出，扩散场法计算HVSR曲线的实质是通过格林函数虚部表示扩散场中包含散射波的能量，计算扩散场中一点的水平与竖向方向能量密度之比。因此，扩散场法不再局限于层状半空间的假设，可以考虑场地中散射体对HVSR曲线的影响。根据式（5）和式（6）可知，局部场地上的HVSR曲线可以通过两种方式得到：① 在存在局部地形的二维半空间表面施加水平方向和竖向的单位简谐荷载，计算相应的广义格林函数虚部，进而通过格林函数虚部与能量密度的正比关系得到；② 根据扩散场理论不同波形具有的能量均分性质，通过大量入射固定能量密度比的各个波形的波构造扩散场，计算地表位移互相关得到。这两种计算HVSR曲线的等价性已经由Perton等（2009）在计算半空间的能量波动情况时证明。相较于通过设置大量入射波构造散射场并计算其散射的第二种方式，直接计算广义格林函数虚部的第一种方式计算量较小。本文采用第一种方式对局部地形上的HVSR曲线进行模拟，计算模型如图1所示，D1为基岩上层状半空间区域，D2为层状半空间中内嵌的沉积区域。在局部地形表面施加竖向荷载与水平荷载，以计算对应的格林函数虚部。
2. 采用IBEM方法模拟HVSR曲线
为了求解局部地形在地表单位简谐荷载作用下产生的格林函数虚部，采用波场分离的思想将总波场分解为不考虑局部地形的自由场和由于局部地形存在而产生的散射场。自由场即不存在局部地形的层状半空间中由外荷载产生的波场，利用直接刚度法即可求得。散射场即入射波在局部地形边界上发生散射产生的波场，通过将边界离散为单元并采用斜线格林函数进行模拟。将求得的自由场和散射场位移、牵引力结果代入边界单元的位移与牵引力平衡方程，即可最终确定总波场。
首先对自由场进行求解. 通过求解平面内波动方程，可得到表示层状半空间外部荷载和层间位移关系的刚度矩阵
$$ {\left[ {{{\tilde {\bar R}}_{x1}}{\text{，}}{\text{i}}{{\tilde {\bar R}}_{z1}}{\text{，}} \cdots {\text{，}}{{\tilde {\bar R}}_{xi}}{\text{，}}{\text{i}}{{\tilde {\bar R}}_{zi}}{\text{，}} \cdots {\text{，}}{{\tilde {\bar R}}_{x{\text{（}}n {\text{＋}} 1{\text{）}}}}{\text{，}}{\text{i}}{{\tilde {\bar R}}_{z{\text{（}}n {\text{＋}} 1{\text{）}}}}} \right]^{\rm{T}}} {\text{＝}} {{{\boldsymbol{S}}{{_{{\rm{P}}{\text{，}}{\rm{SV}}}}}}} {\left[ {{{\tilde {\bar u}}_{x1}}{\text{，}}{\text{i}}{{\tilde {\bar u}}_{z1}}{\text{，}} \cdots {\text{，}}{{\tilde {\bar u}}_{xi}}{\text{，}}{\text{i}}{{\tilde {\bar u}}_{zi}}{\text{，}} \cdots {\text{，}}{{\tilde {\bar u}}_{x{\text{（}}n {\text{＋}} 1{\text{）}}}}{\text{，}}{\text{i}}{{\tilde {\bar u}}_{z{\text{（}}n {\text{＋}} 1{\text{）}}}}} \right]^{\rm{T}}}{\text{，}} $$ (8) 式中上标“~”和“−”分别表示对时间t和空间分量x的傅里叶变换；
${\tilde{ \bar R}_{xi}}$ 和${\tilde {\bar R}_{zi}}$ 分别表示施加在第“i”个土层界面上的水平和竖向外荷载幅值；$ {\tilde u_{xi}} $ 和$ {\tilde u_{zi}} $ 分别表示施加在第“i”个土层界面上的水平和竖向位移幅值；SP， SV为层状半空间整体动力刚度矩阵，具体元素可参见Wolf （1985）。在计算集中荷载产生的位移响应时，荷载施加点处的格林函数实部是奇异的，虚部是收敛的。但在采用IBEM方法求解局部地形的HVSR曲线时，需要得到在边界单元上实部与虚部均收敛的位移解。为了避免施加集中荷载造成边界单元上的位移奇异性，在半空间表面施加作用宽度为2
$ \Delta b $ 的等效均布荷载F1（x）和F3（x）代替集中荷载，二者在$ \Delta b $ 足够小的情况下等效。在二维平面内直角坐标系下，将均布荷载F1（x）和F3（x）转换到波数域，写为关于exp（ikx）的傅里叶积分，即$$ {F_j}{\text{（}}k{\text{）}} {\text{＝}} \frac{{{F_0}}}{{2\pi }}\int_{ - \Delta b}^{{} \Delta b} {\exp {}{\text{（}}{\rm{i}}kx{\text{）}}{}{\rm{d}}x} {\text{＝}} \frac{{{F_0}}}{{\pi k}}\sin k\Delta b{\text{，}}\begin{array}{*{20}{c}} {}&{j {\text{＝}} 1{\text{，}}3{\text{，}}} \end{array} $$ (9) 式中，j取1或3时表示水平和垂直方向对地表施加单位荷载情况，F0为简谐荷载幅值，半空间表面上的水平外荷载幅值
${\tilde {\bar R}_{x1}} {\text{＝}} {\text{（}}{\pi k} {\text{）}}^{-1}\sin k\Delta b$ ，其余外荷载幅值为零，外荷载向量表示为${\left[ {{{\tilde {\bar R}}_{x1}},0, \cdots ,0,0} \right]^{\rm{T}}}$ 。同样，对于地表施加竖向单位荷载情况，${\tilde {\bar R}_{z1}} {\text{＝}}{\text{（}}{\pi k}{\text{）}}^{-1}\sin k\Delta b$ ，外荷载向量可表示为${\left[ {0,{{\tilde {\bar R}}_{z1}}, \cdots ,0,0} \right]^{\rm{T}}}$ 。将外荷载向量代入式（8）即可得到土层任意位置的自由场水平和竖向位移${\tilde {\bar u}_{x1}}$ 和${\tilde {\bar u}_{z1}}$ 。以上刚度矩阵的求解是在波数域内进行的，最终还需将求得的位移进行逆傅里叶变换转换到空间域，即$$ \tilde u{\text{（}}x{\text{，}}\omega {\text{）}} {\text{＝}} \int_{ - \infty }^{ {\text{＋}} \infty } {\tilde {\bar u}{\text{（}}k{\text{，}}\omega {\text{）}}} {{\rm{exp}} {{\text{（}}\text{i}}\omega t{\text{）}}}{\rm{d}}k {\text{.}}$$ (10) 对于三维情况，集中荷载Fj可以表示为
$$ {F_j} {\text{＝}} \delta {\text{（}}x{\text{）}}\delta {\text{（}}y{\text{）}}\exp {\text{（}}{\rm{i}}\omega t{\text{）}}{\text{，}}\begin{array}{*{20}{c}} {}&{{}j {\text{＝}} 1{\text{，}}2{\text{，}}3{\text{，}}} \end{array} $$ (11) 式中，δ为狄拉克函数。采用柱坐标系并展开到空间波数域中，荷载Fj可以表示为
$$ {F_j} {\text{＝}} \frac{{{F_0}}}{{4{\pi ^2}}}\int_0^{2\pi } {\int_0^{ {\text{＋}} \infty } {\exp \left[ {{\rm{i}}kr\cos {\text{（}}\theta {\text{－}} \varphi {\text{）}}} \right]k{\rm{d}}\varphi {\rm{d}}r} } {\text{，}} \,\,\,\,\,{}j{\text{＝}} 1{\text{，}}2{\text{，}}3{\text{，}} $$ (12) 式中，φ和r分别为柱坐标系的坐标变量。
与二维情况相同，集中荷载可以等效为作用在半径为
$ \Delta a $ 的小圆盘上的均布荷载，即$$ {F_j}{\text{（}}k{\text{）}} {\text{＝}} - \frac{{{F_0}\Delta a}}{k}{{\rm{J}}_1}{\text{（}}k\Delta a{\text{）}}{\text{，}}\,\,\,\,\,{}j {\text{＝}} 1{\text{，}}2{\text{，}}3{\text{，}} $$ (13) 式中J1表示第一类一阶贝塞尔函数。将荷载代入三维刚度矩阵求解，可得到层状半空间的自由场位移响应。利用自由场地表水平与竖向位移结果，可计算三维层状半空间的HVSR曲线。
为了求解局部地形产生的散射场，采用平面内斜线格林函数对其进行模拟。该函数最先由Wolf （1985）提出，即作用在土层内部的单位斜线荷载对土层任意位置产生的动力响应。散射场的位移和应力可以由格林函数表示为
$$ \left\{ {\begin{array}{*{20}{c}} {u_x^{\rm{g}}{\text{（}} x {\text{）}}} \\ {u_z^{\rm{g}}{\text{（}} x {\text{）}}} \end{array}} \right\} {\text{＝}} {{{\boldsymbol{g}}_u}} {{{\boldsymbol{P}}}} {\text{＝}} \left\{ {\begin{array}{*{20}{c}} {\displaystyle\sum\limits_{i {\text{＝}} 1}^K {{\left[ {{\boldsymbol{g}}_{u_{11}}{\text{（}} {x{\text{，}}{\xi _i}} {\text{）}}p{\text{（}}{\xi _i}{\text{）}} {\text{＋}} {\boldsymbol{g}}_{u_{13}}{\text{（}} {x{\text{，}}{\xi _i}} {\text{）}}r{\text{（}}{\xi _i}{\text{）}}} \right]} {}} } \\ {\displaystyle\sum\limits_{i {\text{＝}} 1}^K {{\left[ {{\boldsymbol{g}}_{u_{31}}{\text{（}} {x{\text{，}}{\xi _i}} {\text{）}}p{\text{（}}{\xi _i}{\text{）}} {\text{＋}} {\boldsymbol{g}}_{u_{33}}{\text{（}} {x{\text{，}}{\xi _i}} {\text{）}}r{\text{（}}{\xi _i}{\text{）}}} \right]} {}} } \end{array}} \right\} {\text{，}}$$ (14) $$ \left\{ {\begin{array}{*{20}{c}} {\sigma_x^{\rm{g}}{\text{（}} x {\text{）}}} \\ {\sigma_z^{\rm{g}}{\text{（}} x {\text{）}}} \end{array}} \right\} {\text{＝}} {{{\boldsymbol{g}}_{\sigma} }} {{{\boldsymbol{P}}}} {\text{＝}} \left\{ {\begin{array}{*{20}{c}} {\displaystyle\sum\limits_{i {\text{＝}} 1}^K {{\left[ {{\boldsymbol{g}}_{\sigma_{11}}{\text{（}} {x{\text{，}}{\xi _i}} {\text{）}}p{\text{（}}{\xi _i}{\text{）}} {\text{＋}} {\boldsymbol{g}}_{\sigma_{13}}{\text{（}} {x{\text{，}}{\xi _i}} {\text{）}}r{\text{（}}{\xi _i}{\text{）}}} \right]} {}} } \\ {\displaystyle\sum\limits_{i {\text{＝}} 1}^K {{\left[ {{\boldsymbol{g}}_{\sigma_{31}}{\text{（}} {x{\text{，}}{\xi _i}} {\text{）}}p{\text{（}}{\xi _i}{\text{）}} {\text{＋}} {\boldsymbol{g}}_{\sigma_{33}}{\text{（}} {x{\text{，}}{\xi _i}} {\text{）}}r{\text{（}}{\xi _i}{\text{）}}} \right]} {}} } \end{array}} \right\} {\text{，}}$$ (15) 式中：
${{\boldsymbol{g}}_{{u_{ij}}}}{\text{（}}x,{\xi _i}{\text{）}}$ 和${{\boldsymbol{g}}_{{\sigma _{ij}}}}{\text{（}}x,{\xi _i}{\text{）}}$ 分别为位移和应力斜线格林函数，i，j取1或3表示水平或垂直方向；K为边界上划分出的边界单元的个数；P为虚拟荷载密度向量；p（ξi）和r（ξi）分别为水平方向和竖直方向的虚拟荷载密度。为了求解局部地形边界上的虚拟荷载密度，将沉积内外边界单元上的位移和牵引力代入平衡方程。当计算点位于沉积范围内时，荷载施加在沉积表面，沉积内部区域D2波场包含自由场和散射场，沉积外部区域D1仅存在散射场。在沉积边界上，沉积内外位移牵引力连续，边界条件可以写为$$ \begin{array}{*{20}{l}} {\left\{ {\begin{array}{*{20}{l}} {\displaystyle\int_s {\left[ {{t^{\rm{g}}_{2x}}{\text{（}} s {\text{）}}{\text{＋}}t_{2x}^{\rm{f}}{\text{（}} s {\text{）}}} \right]} {\rm{d}}s{\text{＝}}\displaystyle\int_s {{t^{\rm{g}}_{1x}}{\text{（}} s {\text{）}}{\rm{d}}s} } {\text{，}}\\ {\displaystyle\int_s {\left[ {{t^{\rm{g}}_{2z}}{\text{（}} s {\text{）}}{\text{＋}}t_{2z}^{\rm{f}}{\text{（}} s {\text{）}}} \right]} {\rm{d}}s{\text{＝}}\displaystyle\int_s {{t^{\rm{g}}_{1z}}{\text{（}} s {\text{）}}{\rm{d}}s} } {\text{，}}\\ {\displaystyle\int_s {\left[ {{u^{\rm{g}}_{2x}}{\text{（}} s {\text{）}}{\text{＋}}u_{2x}^{\rm{f}}{\text{（}} s {\text{）}}} \right]} {\rm{d}}s{\text{＝}}\displaystyle\int_s {{u^{\rm{g}}_{1x}}{\text{（}} s {\text{）}}{\rm{d}}s} } {\text{，}}\\ {\displaystyle\int_s {\left[ {{u^{\rm{g}}_{2z}}{\text{（}} s {\text{）}}{\text{＋}}u_{2z}^{\rm{f}}{\text{（}} s {\text{）}}} \right]} {\rm{d}}s{\text{＝}}\displaystyle\int_s {{u^{\rm{g}}_{1z}}{\text{（}} s {\text{）}}{\rm{d}}s} }{\text{，}} \end{array}} \right.}&{i {\text{＝}} 1{\text{，}}2{\text{，}} \cdots {\text{，}}K}{\text{，}} \end{array} $$ (16) 式中，t为边界单元上由应力合成得到的牵引力，下标1，2分别代表沉积地形外部和内部，上标g，f分别代表散射场和自由场。当计算点位于沉积范围外时，荷载施加在沉积外部，沉积内部区域D2波场仅包括散射场，沉积外部区域D1波场包括自由场和散射场，边界条件可以写为
$$ \begin{array}{*{20}{l}} {\left\{ {\begin{array}{*{20}{l}} {\displaystyle\int_s {\left[ {{t^{\rm{g}}_{1x}}{\text{（}} s {\text{）}}{\text{+}}t_{1x}^{\rm{f}}{\text{（}} s {\text{）}}} \right]} {\rm{d}}s{\text{＝}}\displaystyle\int_s {{t^{\rm{g}}_{2x}}{\text{（}} s {\text{）}}{\rm{d}}s} }{\text{，}} \\ {\displaystyle\int_s {\left[ {{t^{\rm{g}}_{1z}}{\text{（}} s {\text{）}}{\text{+}}t_{1z}^{\rm{f}}{\text{（}} s {\text{）}}} \right]} {\rm{d}}s{\text{＝}}\displaystyle\int_s {{t^{\rm{g}}_{2z}}{\text{（}} s {\text{）}}{\rm{d}}s} }{\text{，}} \\ {\displaystyle\int_s {\left[ {{u^{\rm{g}}_{1x}}{\text{（}} s {\text{）}}{\text{+}}u_{1x}^{\rm{f}}{\text{（}} s {\text{）}}} \right]} {\rm{d}}s{\text{＝}}\displaystyle\int_s {u_{2x}^{{\rm{g}}}{\text{（}} s {\text{）}}{\rm{d}}s} }{\text{，}} \\ {\displaystyle\int_s {\left[ {{u^{\rm{g}}_{1z}}{\text{（}} s {\text{）}}{\text{+}}u_{1z}^{\rm{f}}{\text{（}} s {\text{）}}} \right]} {\rm{d}}s{\text{＝}}\displaystyle\int_s {{u^{\rm{g}}_{2z}}{\text{（}} s {\text{）}}{\rm{d}}s} } {\text{，}} \end{array}} \right.}&{i {\text{＝}} 1{\text{，}}2{\text{，}} \cdots {\text{，}}K} {\text{.}} \end{array} $$ (17) 求解式（16）和（17）可得局部地形边界单元上的荷载密度，确定散射波场。将得到的自由场与散射场叠加，即可得到总波场。对地表计算点的总场结果（广义格林函数）取虚部，即可得到式（7）中局部地形场地上的
${{\rm{Im}}} \left[ {{G_{11}}{\text{（}}x{\text{，}}x{\text{，}}\omega {\text{）}}} \right]$ 和${{\rm{Im}}} \left[ {{G_{33}}{\text{（}}x{\text{，}}x{\text{，}}\omega {\text{）}}} \right]$ ，进而得到计算点处的HVSR曲线。3. 方法验证
由于半空间中散射体的不均匀分布，波场可能不具有扩散场性质（Perton et al，2009）。因此，需要对二维层状半空间中沉积地形的扩散场理论的适用性进行验证。Perton等（2009）指出，波场的扩散场性质与能量均分性质等价，因此对场地的能量均分进行验证即可证明波场为扩散场。
假设层状半空间中沉积地形场地符合扩散场能量均分，采用式（3）对其沿深度的能量波动情况进行计算。将局部地形场地取为层状半空间中梯形沉积地形，梯形沉积下底宽为40 m，上底宽为80 m，沉积深度为40 m，基岩上单一土层厚度取为80 m。沉积、外部土层和基岩材料阻抗比为1∶2∶4，泊松比均取为0.25，阻尼比取为0.005。则沉积地形中心点处的方向能量密度随深度变化如图2所示.可以看出，沉积中心点位置处，在沉积内部范围内（0—40 m）和沉积外部土层范围内（40—80 m），水平方向和竖直方向的能量密度在同一数值范围内波动。随着深度增加，水平方向和竖直方向的能量密度逐渐趋于相等（图2a）；在梯形沉积地形的非沉积中心点，竖直、水平两方向的能量密度基本在同一数值范围内波动（图2b）。但对比中心点处的能量波动情况，可以看出在阻抗比界面（深度为40 m）附近，非中心点处的水平与竖向能量密度发生了显著的波动。该现象是由波在横向的阻抗比界面上发生反射造成的。由于受到横向阻抗比界面的影响范围较小，沉积地形非中心点的能量水平总体仍符合扩散场能量均分特征，即波场的各方向能量密度相等（Sánchez-Sesma et al，2008），满足使用式（7）计算HVSR曲线的条件。
![]()
图 2 层状半空间中沉积内部中心点x＝0 m （a） 和非中心点x＝16 m （b）的方向能量密度随深度变化情况Figure 2. Energy fluctuation of center point x＝0 m （a） and non-central point x＝16 m （b）in depth direction inside the alluvial canyon in layered half-space其次，为了验证采用本文方法得到结果的正确性，与已有结果进行对比验证。首先采用直接刚度法，计算三维情况下层状半空间HVSR曲线并与Sánchez-Sesma等（2011）给出的结果进行对比，验证自由场与层状半空间求解结果的正确性。层状半空间土层厚度40 m，具体参数列于表1。
表 1 层状半空间计算参数Table 1. Calculation parameters of layered half space剪切波速
vS/（m·s−1）泊松比ν 土层密度
ρ/（kg·m−3）阻尼比ζ 土层 70 0.496 1 200 0.05 基岩 1 000 0.333 2 500 0.05 本文采用刚度矩阵方法计算得到的HVSR曲线与Sánchez-Sesma等（2011）结果吻合良好（图3），可以证明自由场与层状半空间HVSR曲线求解的正确性。
![]()
图 3 采用本文方法得到的层状半空间HVSR曲线与Sánchez-Sesma等（2011）结果对比Figure 3. Comparisons of the result of layered half-space HVSR curves in this method with Sánchez-Sesma et al （2011）目前尚未有研究给出采用扩散场理论计算的二维平面内局部地形场地的HVSR曲线结果，仅Perton和Sánchez-Sesma （2016）通过扩散场理论计算了均匀半空间中沉积地形表面的格林函数张量。为了验证采用IBEM 方法得到的二维沉积地形上格林函数虚部的正确性，将本文结果与Perton和Sánchez-Sesma （2016）结果进行对比， 结果如图4所示。可以看出两结果总体对应良好，本文方法的正确性得以验证。
![]()
图 4 二维平面内采用本文方法得到的沉积地形上格林函数张量与Perton和Sánchez-Sesma （2016）的结果对比Figure 4. Comparisons of the result of Green's function tensor on sedimentary topography in this method with Perton and Sánchez-Sesma （2016）4. 算例计算与结果分析
本节首先通过对比层状半空间与二维沉积地形的HVSR曲线，分析HVSR曲线对层状半空间参数变化和沉积地形存在的敏感性；其次，通过对沉积地形的不同特征，包括沉积内外材料阻抗比、沉积形状（梯形沉积侧面边界角度）以及计算点与沉积的相对位置对HVSR曲线的影响规律进行分析；最后，基于算例分析得出的HVSR曲线特点，探讨通过HVSR曲线判断沉积地形的大小、沉积材料情况和相对位置的可能性。本节算例计算中，将沉积地形形状取为梯形，计算模型如图5所示。
![]()
图 5 梯形沉积地形计算模型图Figure 5. Calculation model graph of trape-zoidal sedimentary topography4.1 层状半空间与二维沉积地形HVSR曲线对比
层状半空间的HVSR曲线的第一峰值频率与土层的卓越频率一致。但在对比分析层状半空间与沉积地形的HVSR曲线之前，仍需明确层状半空间速度结构对HVSR曲线的幅值和形态的影响规律。首先，分别计算四个不同沉积土层厚度、阻抗比和土层数的层状半空间HVSR曲线，层状半空间的土层厚度和材料类别如图6所示。不同沉积土层与基岩材料参数列于表2，得出相应层状半空间的HVSR曲线计算结果（图7）。
![]()
图 6 层状半空间①—④ （a−d）计算模型图Figure 6. Calculation model graph of layered half-spaces ①—④ （a−d）表 2 层状半空间不同阻抗比情况计算参数Table 2. Parameters of different impedance ratios in layered half-space剪切波速
vS/（m·s−1）泊松比
ν土层密度
ρ/（kg·m−3）阻尼比
ζ沉积土层① 140 0.496 1 200 0.05 沉积土层② 280 0.496 1 200 0.05 基岩半空间① 560 0.496 1 200 0.05 基岩半空间② 280 0.496 1 200 0.05 从图7中可以看出，不同层状半空间因其具有不同的阻抗比和土层厚度，故对HVSR曲线具有较大影响。对比层状半空间③（图6c）和②（图6b）可以看出，相同土层厚度情况下，较高的阻抗比界面会使HVSR曲线峰值点幅值增大。层状半空间①（图6a）两个阻抗比界面）对比其它层状半空间（单一阻抗比界面）的HVSR曲线可以看出，多个阻抗比界面会引起第一峰值点幅值的变化，但并没有改变HVSR曲线第一峰值显著、高频段曲线幅值逐渐减小这一基本形态。
为了对比层状半空间和二维沉积地形的HVSR曲线，计算均匀半空间中梯形沉积地形的HVSR曲线，并与计算点处对应的层状半空间结果进行对比，将层状半空间基岩上土层厚度D取为60 m，沉积厚度a取为30 m，沉积内外材料阻抗比为1∶2，结果如图8所示。可以看出，沉积地形上的HVSR曲线与层状半空间结果差异显著，沉积地形两侧的阻抗界面使得HVSR曲线的第一峰值明显削弱，而第二峰值得到了增强。HVSR曲线形态对沉积地形的存在更为敏感，不再符合层状半空间情况下第一峰值点的幅值最大，后续峰值幅值逐渐减小的规律，而是第一峰值附近出现“平台”现象，第二峰值更为突出，此现象与Uebayashi （2003）所给出的结论一致。
![]()
图 8 沉积地形与相应层状半空间的HVSR曲线Figure 8. HVSR curves of sedimentary topography and corresponding layered half-space本节计算结果表明，土层的不同阻抗比仅会使得层状半空间的HVSR曲线峰值点幅值和频率在一定范围内变化，不会改变HVSR曲线的基本形态。而沉积地形的存在则显著地改变了HVSR曲线的特征，在后续分析中，可认为HVSR曲线的形态改变是由沉积地形的存在造成的。
4.2 沉积内外材料阻抗比对HVSR曲线的影响
层状半空间的土层与基岩阻抗比会显著改变层状半空间上HVSR曲线的形态，因此在具有明显阻抗比的沉积地形上，HVSR曲线也会受沉积内外介质阻抗比大小的影响。为了研究其影响规律，将沉积外部土层阻抗取为沉积内部材料的2倍、3倍、4倍和5倍。场地模型尺寸与4.1节中沉积算例相同，沉积地形参数列于表3。
表 3 不同沉积内外材料阻抗比情况计算参数Table 3. Parameters of the alluvial canyon materials with different impedance ratios剪切波速
vS/（m·s−1）泊松比
ν土层密度
ρ/（kg·m−3）阻尼比
ζ沉积内部土层 140 0.496 1 200 0.05 外部土层① 280 0.496 1 200 0.05 外部土层② 420 0.496 1 200 0.05 外部土层③ 560 0.496 1 200 0.05 外部土层④ 700 0.496 1 200 0.05 基岩半空间 1 000 0.333 2 500 0.05 不同阻抗比的局部地形与计算点下对应的一维层状半空间的HVSR曲线如图9所示。对比沉积地形（图9a）与相应层状半空间的HVSR曲线（图9b），可以看出二者存在明显的差异。沉积地形的存在使得HVSR曲线的第一峰值点相对层状半空间情况幅值减小，频率增大，第二峰值点幅值增大，频率增大，曲线形态显著改变。同时，沉积地形对HVSR曲线的影响随着沉积内外阻抗差异的增大而增大。在HVSR曲线幅值方面，第一峰值点的幅值随阻抗差异增大而减小，第二峰值点的幅值随阻抗差异增大而增大。在HVSR曲线峰值频率方面，第一、第二峰值频率均随阻抗差异增大而增大，但第一峰值频率受局部地形影响增大更为明显。对于算例中的沉积地形，当沉积内外阻抗比达到1 ∶ 5时，原有第一峰值逐渐消失，现有第一峰值频率达到相应层状半空间HVSR曲线的3.3倍，峰值点幅值仅为层状半空间幅值的28.3％。同时阻抗比1 ∶ 5情况下的曲线第一峰值频率是阻抗比1 ∶ 2情况下的3.2倍，峰值点幅值仅为阻抗比1 ∶ 2情况下的65.3％。在HVSR曲线形态方面，随着沉积内外阻抗差异的增大，HVSR曲线第一峰值点幅值逐渐减小，形成平台现象，第二峰值点幅值逐渐增大。
![]()
图 9 沉积地形（a）和其对应的层状半空间（b）内不同沉积内外阻抗比的HVSR曲线Figure 9. HVSR curves of the sedimentary materials with different impedance ratios for sedimentarytopography (a) and its corresponding layered half space (b)4.3 沉积形状变化对HVSR曲线的影响
为了探究局部地形的形状（横向不均匀界面坡度）对HVSR曲线的影响，本研究计算了具有不同侧面边界坡度的梯形沉积上的HVSR曲线并进行对比。计算模型中沉积厚度仍取为基岩上土层厚度D＝60 m，沉积厚度a＝30 m。梯形沉积侧面边界与水平方向夹角α分别取15°，30°，45°，60°和75°。沉积内部材料、外部土层材料以及基岩材料参数列于表4，将沉积外部土层材料参数分别取为表4中外部土层①与外部土层②，计算层状半空间中沉积内外阻抗比1∶2和1∶4情况下的HVSR曲线。
表 4 不同沉积形状情况沉积地形计算参数Table 4. Parameters of the alluvial canyon materials with different topography shapes剪切波速
vS/（m·s−1）泊松比
ν土层密度
ρ/（kg·m−3）阻尼比
ζ沉积内部土层 140 0.496 1 200 0.05 外部土层① 280 0.496 1 200 0.05 外部土层② 420 0.496 1 200 0.05 基岩半空间 1 000 0.333 2 500 0.05 如图10所示，沉积形状显著影响了HVSR曲线的峰值点频率和幅值。对比层状半空间情况与α为15°时的HVSR曲线可以看出，在坡度较小的情况下，局部地形的HVSR曲线与层状半空间结果十分接近。随着沉积地形界面坡度增大，曲线第一峰值点频率随沉积地形界面坡度的增大而增大，曲线第一峰值点幅值随沉积地形界面坡度增大而减小。从图10b中可以看出，在沉积地形界面逐渐接近竖直的情况下，HVSR曲线第一峰值点逐渐向右下方移动，最终第一峰值逐渐消失，出现“平台”现象.当沉积内外阻抗比为1∶4，α＝15°时，HVSR曲线第一峰值点幅值是层状半空间的1.54倍，当α＝60°时，HVSR曲线第一峰值点频率是层状半空间的1.77倍。对比阻抗比为1∶2 （图10a）和1∶4 （图10b）的HVSR曲线，在阻抗差异较大情况下，沉积地形的形状对HVSR曲线影响更显著。
![]()
图 10 沉积内外材料阻抗比为1∶2 （a）和1∶4 （b）时不同横向不均匀界面坡度的HVSR曲线Figure 10. HVSRs of different laterally inhomogeneous interface slopes with impedance ratios 1∶2 （a） and 1∶4 （b） between inside and outside materials of sedimental topography4.4 不同计算点相对位置的HVSR曲线
对于同一沉积地形，不同计算点位置的HVSR曲线显著不同。因此，本小节基于4.1节中的计算模型，将单位简谐荷载分别施加在沉积表面x/a为0，0.2，0.4，0.6和0.8的位置处，计算相应位置的HVSR曲线。将沉积内土层材料参数取为与4.1节中沉积内部土层相同，并将沉积外部土层分别取为4.1节中外部土层①和外部土层②材料，计算沉积内外材料阻抗比分别为1∶2和1∶4的HVSR曲线（图11）。
![]()
图 11 沉积表面不同相对位置x/a上的HVSR曲线（a） 沉积内外材料阻抗比1∶2；（b） 沉积内外材料阻抗比1∶4Figure 11. HVSR curves of the sedimentary materials with different positions x/a（a） The impedance ratio of internal and external material is 1∶2；（b） The impedance ratio of internal and external material is 1∶4如图11所示，随着计算点逐渐靠近沉积边界（x/a增大），沉积地形的存在使得HVSR曲线形态与沉积中心位置（x/a＝0.0处）对应的层状半空间情况结果差异增大.当沉积内外材料阻抗比为1∶2时，随着计算点到沉积边界距离减小，HVSR曲线第一峰值点幅值和第二峰值点幅值逐渐减小，峰值频率基本不变。沉积表面的HVSR曲线第一峰值点处幅值在层状半空间HVSR第一峰值的37.4%—49.2%范围内变化（图11a）；在沉积内外材料阻抗比为1∶4时，随着计算点到沉积边界距离减小，HVSR曲线第一峰值点幅值逐渐增大，第二峰值点幅值逐渐减小，峰值点频率基本不变。沉积表面的HVSR曲线第一峰值点处幅值在层状半空间HVSR的第一峰值的29.3%—50.9%范围内变化（图11b）。对比图11a和图11b，可以看出在阻抗差异较大的情况下，沉积地形不同位置的HVSR曲线与层状半空间情况差异更为明显，尤其是距离边界较近位置上的HVSR曲线形态与层状半空间情况明显不同。在距离沉积边界较近的计算点，更容易接收到沉积边界反射的高频波。因此在第二峰值后的中高频段，HVSR曲线幅值相对较大。随着计算点到沉积边界距离的减小，距离沉积地形边界较近位置处的HVSR曲线的形态变为第二峰值附近频段幅值相对较小，而中高频段幅值相对较大。
综合考虑局部地形对HVSR曲线的幅值和形态的影响，距离边界越近，高频段的HVSR曲线波动越剧烈，并且随着频率增大高频段曲线幅值逐渐增大。该结果在阻抗比和沉积形状对HVSR曲线的影响规律中并未显现。
4.5 根据HVSR曲线判断局部地形特征
上述计算初步分析了沉积地形对HVSR曲线的影响规律，结果表明，沉积地形的存在使得HVSR曲线明显区别于层状半空间结果，这使得根据HVSR曲线特征勘探未知场地的局部地形分布情况成为可能。同时，HVSR方法与传统场地勘探方法相比成本较低，测量时间较短，具有明显的成本优势。本节通过总结分析沉积地形上HVSR曲线的规律，对HVSR方法勘探局部地形的可行性进行探讨。
在获得场地实测的HVSR曲线后，需要根据HVSR曲线的峰值频率、幅值以及曲线形态等特征筛选出与层状半空间曲线差异明显的计算点，初步判断局部地形位置。HVSR曲线的第一、第二峰值点的频率和幅值均会受到沉积地形的阻抗差异、沉积形状和计算点位置的影响。虽然单一因素对其具有较明显的影响规律，但在多因素共同作用下，HVSR曲线的第一、第二峰值较难作为判断局部地形情况的依据。而局部场地上的HVSR曲线的形态具有鲜明的特点。若所测得的HVSR曲线存在显著的“平台”现象或高频幅值较大现象，则可以认为测点附近存在局部地形。
在根据HVSR曲线形态确定局部地形的大体位置后，尝试对局部地形的细节特征进行判断。若测点附近已有钻孔资料，则可以得出测点下的土层一维速度结构，计算出一维层状半空间的HVSR曲线。理论上，通过一维层状半空间的HVSR曲线与实测曲线相对比，结合局部地形单一特征对HVSR曲线的影响规律，可以对沉积特征进行判断。但由于局部地形的阻抗比和界面坡度对HVSR曲线的影响是耦合的，因此不易根据HVSR曲线的峰值点频率和幅值将局部地形的特征分辨出来。若已有测点附近的地震剖面资料，则可以得到局部地形界面坡度情况，通过数值计算与实测HVSR曲线的对比来判断局部地形的其它细节。
综上所述，采用HVSR方法能够对层状场地中的局部地形进行探测，以HVSR曲线形态为依据判断测点附近是否存在横向不均匀体，但对局部地形的形状、阻抗比等具体特征仅凭HVSR曲线难以确定。因此，在勘探地下结构未知的场地局部地形分布情况时，可采用HVSR方法和其它勘探方法相结合的方式，首先采用成本较低的HVSR方法进行初次勘探，确定局部地形位置，再通过其它勘探手段对已确定存在局部地形的区域进行细部勘探。引入HVSR方法进行联合勘探可以节约大量勘探时间和成本。
5. 讨论与结论
HVSR曲线可以较为准确地反映土层的卓越频率，具有较好的研究意义和应用价值。但大多数研究未考虑局部地形的影响，忽视了局部地形对HVSR曲线的研究和应用造成的干扰。本文采用Kawase等（2011）和Sánchez-Sesma等（2011）提出的扩散场方法，结合IBEM方法计算了沉积地形中的HVSR曲线，对沉积内外材料阻抗比、沉积形状以及计算点位置对HVSR曲线的影响进行了分析。结果表明，沉积地形的存在会显著影响附近位置处的HVSR曲线形态，具体体现在以下几个方面：
1） HVSR曲线对沉积地形内外的材料阻抗比变化最为敏感，沉积地形形状和计算点位置对HVSR曲线的影响均与阻抗比有关，沉积内外阻抗差异越大，沉积地形对HVSR曲线的影响越显著；
2） 在沉积内外阻抗差异较大的情况下，沉积地形的存在使得HVSR曲线的第一峰值频率相对于层状半空间情况显著增大，第一峰值点幅值显著减小，HVSR曲线形态出现平台现象；
3） 沉积地形形状对HVSR曲线的影响主要体现为第一峰值点处频率随阻抗界面坡度增大而增大，幅值随阻抗界面坡度增大而减小；
4） 随着与沉积边界距离的减小，不同计算点位置上的HVSR曲线形态明显改变，同时高频段曲线幅值相对于层状半空间情况增大。
在上述基础上，本文探讨了采用HVSR方法勘探场地局部地形特征的可行性。局部地形上的HVSR曲线特征表明，可采用HVSR方法进行初步勘探确定场地局部地形分布位置，再结合其它勘探方法确定局部地形的细部情况，从而大大节约局部地形场地的勘探成本。
需要指出的是，本文尚未解决局部地形特征对HVSR曲线的耦合影响问题，无法通过HVSR曲线判断局部地形的细部特征，仅就局部地形之一的各向同性沉积地形进行了HVSR曲线计算与研究。在后续研究中，本文方法还可以应用于计算其它形式的2.5维或三维局部地形以及考虑场地土层各向异性的情况。若能采用本文方法模拟出实际地形的HVSR曲线，则有可能通过大量的数值计算反演出场地局部地形的细部特征，发挥HVSR方法潜在的应用价值。
-
![]()
图 2 层状半空间中沉积内部中心点x＝0 m （a） 和非中心点x＝16 m （b）的方向能量密度随深度变化情况
Figure 2. Energy fluctuation of center point x＝0 m （a） and non-central point x＝16 m （b）in depth direction inside the alluvial canyon in layered half-space
![]()
图 3 采用本文方法得到的层状半空间HVSR曲线与Sánchez-Sesma等（2011）结果对比
Figure 3. Comparisons of the result of layered half-space HVSR curves in this method with Sánchez-Sesma et al （2011）
![]()
图 4 二维平面内采用本文方法得到的沉积地形上格林函数张量与Perton和Sánchez-Sesma （2016）的结果对比
Figure 4. Comparisons of the result of Green's function tensor on sedimentary topography in this method with Perton and Sánchez-Sesma （2016）
![]()
图 5 梯形沉积地形计算模型图
Figure 5. Calculation model graph of trape-zoidal sedimentary topography
![]()
图 6 层状半空间①—④ （a−d）计算模型图
Figure 6. Calculation model graph of layered half-spaces ①—④ （a−d）
![]()
图 8 沉积地形与相应层状半空间的HVSR曲线
Figure 8. HVSR curves of sedimentary topography and corresponding layered half-space
![]()
图 9 沉积地形（a）和其对应的层状半空间（b）内不同沉积内外阻抗比的HVSR曲线
Figure 9. HVSR curves of the sedimentary materials with different impedance ratios for sedimentarytopography (a) and its corresponding layered half space (b)
![]()
图 10 沉积内外材料阻抗比为1∶2 （a）和1∶4 （b）时不同横向不均匀界面坡度的HVSR曲线
Figure 10. HVSRs of different laterally inhomogeneous interface slopes with impedance ratios 1∶2 （a） and 1∶4 （b） between inside and outside materials of sedimental topography
![]()
图 11 沉积表面不同相对位置x/a上的HVSR曲线
（a） 沉积内外材料阻抗比1∶2；（b） 沉积内外材料阻抗比1∶4
Figure 11. HVSR curves of the sedimentary materials with different positions x/a
（a） The impedance ratio of internal and external material is 1∶2；（b） The impedance ratio of internal and external material is 1∶4
表 1 层状半空间计算参数
Table 1 Calculation parameters of layered half space
剪切波速
vS/（m·s−1）泊松比ν 土层密度
ρ/（kg·m−3）阻尼比ζ 土层 70 0.496 1 200 0.05 基岩 1 000 0.333 2 500 0.05 
下载: 导出CSV
表 2 层状半空间不同阻抗比情况计算参数
Table 2 Parameters of different impedance ratios in layered half-space
剪切波速
vS/（m·s−1）泊松比
ν土层密度
ρ/（kg·m−3）阻尼比
ζ沉积土层① 140 0.496 1 200 0.05 沉积土层② 280 0.496 1 200 0.05 基岩半空间① 560 0.496 1 200 0.05 基岩半空间② 280 0.496 1 200 0.05
下载: 导出CSV
表 3 不同沉积内外材料阻抗比情况计算参数
Table 3 Parameters of the alluvial canyon materials with different impedance ratios
剪切波速
vS/（m·s−1）泊松比
ν土层密度
ρ/（kg·m−3）阻尼比
ζ沉积内部土层 140 0.496 1 200 0.05 外部土层① 280 0.496 1 200 0.05 外部土层② 420 0.496 1 200 0.05 外部土层③ 560 0.496 1 200 0.05 外部土层④ 700 0.496 1 200 0.05 基岩半空间 1 000 0.333 2 500 0.05
下载: 导出CSV
表 4 不同沉积形状情况沉积地形计算参数
Table 4 Parameters of the alluvial canyon materials with different topography shapes
剪切波速
vS/（m·s−1）泊松比
ν土层密度
ρ/（kg·m−3）阻尼比
ζ沉积内部土层 140 0.496 1 200 0.05 外部土层① 280 0.496 1 200 0.05 外部土层② 420 0.496 1 200 0.05 基岩半空间 1 000 0.333 2 500 0.05
下载: 导出CSV
-
林国良,张潜,崔建文,赵昆,杨黎薇. 2019. 利用地脉动HVSR研究2014年鲁甸6.5级地震场地效应[J]. 地震研究,42(4):531–537. doi: 10.3969/j.issn.1000-0666.2019.04.011 Lin G L,Zhang Q,Cui J W,Zhao K,Yang L W. 2019. Determining the site effects of the 2014 Ludian MS6.5 earthquake using HVSR microtremor method[J]. Journal of Seismological Research,42(4):531–537 (in Chinese).
卢育霞,刘琨,王良,魏来,李少华. 2017. 基于台阵记录的土层山体场地效应分析[J]. 地震学报,39(6):941–954. Lu Y X,Liu K,Wang L,Wei L,Li S H. 2017. Site effect of unconsolidated soil hill based on seismic array records[J]. Acta Seismologica Sinica,39(6):941–954 (in Chinese).
荣棉水,李小军,王振明,吕悦军. 2016. HVSR方法用于地震作用下场地效应分析的适用性研究[J]. 地球物理学报,59(8):2878–2891. doi: 10.6038/cjg20160814 Rong M S,Li X J,Wang Z M,Lü Y J. 2016. Applicability of HVSR in an analysis of site-effects caused by earthquakes[J]. Chinese Journal of Geophysics,59(8):2878–2891 (in Chinese).
荣棉水,符力耘,李小军. 2018. 基于单台加速度记录的混合全局优化HVSR反演场地浅层速度结构[J]. 地球物理学报,61(3):938–947. doi: 10.6038/cjg2018L0171 Rong M S,Fu L Y,Li X J. 2018. Inversion of site velocity structure using a hybrid global optimization algorithm based on HVSRs of accelerograms recorded by a single station[J]. Chinese Journal of Geophysics,61(3):938–947 (in Chinese).
王伟君,陈棋福,齐诚,谭毅培,张项,周青云. 2011. 利用噪声HVSR方法探测近地表结构的可能性和局限性:以保定地区为例[J]. 地球物理学报,54(7):1783–1797. doi: 10.3969/j.issn.0001-5733.2011.07.012 Wang W J,Chen Q F,Qi C,Tan Y P,Zhang X,Zhou Q Y. 2011. The feasibilities and limitations to explore the near-surface structure with microtremor HVSR method:A case in Baoding area of Hebei Province,China[J]. Chinese Journal of Geophysics,54(7):1783–1797 (in Chinese).
温瑞智,冀昆,任叶飞,王宏伟. 2015. 基于谱比法的我国强震台站场地分类[J]. 岩石力学与工程学报,34(6):1236–1241. Wen R Z,Ji K,Ren Y F,Wang H W. 2015. Site classification for strong earthquake stations in China using spectral ratio method[J]. Chinese Journal of Rock Mechanics and Engineering,34(6):1236–1241 (in Chinese).
Arai H,Tokimatsu K. 2004. S-wave velocity profiling by inversion of microtremor H/V spectrum[J]. Bull Seismol Soc Am,94(1):53–63. doi: 10.1785/0120030028
Chavez-Gárcía F J,Rodríguez M,Field E,Hatzfeld D. 1997. Topographic site effects:A comparison of two nonreference methods[J]. Bull Seismol Soc Am,87(6):1667–1673. doi: 10.1785/BSSA0870061667
Gosar A. 2007. Microtremor HVSR study for assessing site effects in the Bovec basin (NW Slovenia) related to 1998 MW5.6 and 2004 MW5.2 earthquakes[J]. Eng Geol,91(2/3/4):178–193.
Guéguen P,Cornou C,Garambois S,Banton J. 2007. On the limitation of the H/V spectral ratio using seismic noise as an exploration tool:Application to the Grenoble Valley (France),a small apex ratio basin[J]. Pure Appl Geophys,164(1):115–134. doi: 10.1007/s00024-006-0151-x
Kawase H,Sánchez-Sesma F J,Matsushima S. 2011. The optimal use of horizontal-to-vertical spectral ratios of earthquake motions for velocity inversions based on diffuse-field theory for plane waves[J]. Bull Seismol Soc Am,101(5):2001–2014. doi: 10.1785/0120100263
Matsushima S,Hirokawa T,De Martin F,Kawase H,Sánchez-Sesma F J. 2014. The effect of lateral heterogeneity on horizontal-to-vertical spectral ratio of microtremors inferred from observation and synthetics[J]. Bull Seismol Soc Am,104(1):381–393. doi: 10.1785/0120120321
Nakamura Y. 1989. A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface[J]. Quarterly Reports Railway Tech Res Inst,30(1):25–33.
Napolitano F,Gvasi A,La Rocca M L,Guerra I,Scarpa R. 2018. Site effects in the Pollino region from the HVSR and polarization of seismic noise and earthquakes[J]. Bull Seismol Soc Am,108(1):309–321. doi: 10.1785/0120170197
Nogoshi M,Igarashi T. 1971. On the amplitude characteristics of microtremor:Part 2[J]. Seism Soc Jap,24:26–40.
Perton M,Sánchez-Sesma F J,Rodríguez-Castellanos A,Campillo M,Weaver R L. 2009. Two perspectives on equipartition in diffuse elastic fields in three dimensions[J]. J Acoust Soc Am,126(3):1125–1130. doi: 10.1121/1.3177262
Perton M,Sánchez-Sesma F J. 2016. Green's function calculation from equipartition theorem[J]. J Acoust Soc Am,140(2):1309–1318. doi: 10.1121/1.4961208
Perton M,Spica Z,Caudron C. 2018. Inversion of the horizontal-to-vertical spectral ratio in presence of strong lateral heterogeneity[J]. Geophys J Int,212(2):930–941. doi: 10.1093/gji/ggx458
Piña-Flores J,Perton M,García-Jerez A,Carmona E,Luzón F,Molina-Villegas J C,Sánchez-Sesma F J. 2017. The inversion of spectral ratio H/V in a layered system using the diffuse field assumption (DFA)[J]. Geophys J Int,208(1):577–588. doi: 10.1093/gji/ggw416
Sánchez-Sesma F J,Campillo M. 2006. Retrieval of the Green’s function from cross correlation:The canonical elastic problem[J]. Bull Seismol Soc Am,96(3):1182–1191. doi: 10.1785/0120050181
Sánchez-Sesma F J,Pérez-Ruiz J A,Luzón F,Campillo M,Rodríguez-Castellanos A. 2008. Diffuse fields in dynamic elasticity[J]. Wave Motion,45(5):641–654. doi: 10.1016/j.wavemoti.2007.07.005
Sánchez-Sesma F J,Rodríguez M,Iturrarán-Viveros U,Luzón F,Campillo M,Margerin L,García-Jerez A,Suarez M,Santoyo M A,Rodríguez-Castellanos A. 2011. A theory for microtremor H/V spectral ratio:Application for a layered medium[J]. Geophys J Int,186(1):221–225. doi: 10.1111/j.1365-246X.2011.05064.x
Sánchez-Sesma F J. 2017. Modeling and inversion of the microtremor H/V spectral ratio:Physical basis behind the diffuse field approach[J]. Earth Planets Space,69(1):92. doi: 10.1186/s40623-017-0667-6
Stanko D,Markušić S,Strelec S,Gazdek M. 2017. HVSR analysis of seismic site effects and soil-structure resonance in Varaždin city (North Croatia)[J]. Soil Dyn Earthq Eng,92:666–677. doi: 10.1016/j.soildyn.2016.10.022
Uebayashi H. 2003. Extrapolation of irregular subsurface structures using the horizontal-to-vertical spectral ratio of long-period microtremors[J]. Bull Seismol Soc Am,93(2):570–582. doi: 10.1785/0120020137
Uebayashi H,Kawabe H,Kamae K. 2012. Reproduction of microseism H/V spectral features using a three-dimensional complex topographical model of the sediment-bedrock interface in the Osaka sedimentary basin[J]. Geophys J Int,189(2):1060–1074. doi: 10.1111/j.1365-246X.2012.05408.x
Weaver R L. 1982. On diffuse waves in solid media[J]. J Acoust Soc Am,71(6):1608–1609. doi: 10.1121/1.387816
Weaver R L. 1985. Diffuse elastic waves at a free surface[J]. J Acoust Soc Am,78(1):131–136. doi: 10.1121/1.392576
Wolf J P. 1985. Dynamic Soil-Structure Interaction[M]. Englewood Cliffs: Prentice-Hall.
计量
- 文章访问数: 263
- HTML全文浏览量: 202
- PDF下载量: 43
