近场地震动的周期特性、PGV/PGA、持时等参数对减隔震桥梁和非减隔震桥梁的地震反应均有显著的影响^[1-6]，历史地震记录或者人工波虽然可以实现具有与弹性设计反应谱一致的反应谱特性，但不能直接反映地震断层参数，其周期特性、PGV/PGA和时历特性没有考虑潜在震源个性的影响，有可能导致桥梁地震反应计算结果不合理.

格林函数法根据断层破坏过程以及地震波传播过程模拟观测点地震动，自Hartzell^[7]提出经验格林函数法以来，在地震学领域得到迅速发展，在反演历史地震和强地震动模拟中得到广泛应用. 但是经验格林函数需要震源参数、传播路径和场地条件相同的小震记录，当缺乏小震记录时难以实现场地地震动的模拟. Boore^[8]提出用人工随机波替代小震记录的随机格林函数法，为设计地震动模拟提供了一种途径. 近年来，随机格林函数法模拟设计地震动的方法受到学者们的关注，如陶夏新等^[9]分析近断层强地震动的破裂方向性效应和上盘效应；王国新等^[10]预测沈阳地区近断层地震动的分布特性；王海云等^[11]用3个台站的地震记录校正随机格林函数，模拟 1997 年新疆伽师地震周围的近断层地震动分布；孙吉泽等^[12]模拟2013年乌鲁木齐 M_S5.6 和 M_S5.1 地震的强地震动分布；李启成等^[13]用基于动力学拐角模拟的随机有限断层法模拟2016年新西兰M_W7.8地震的地震动，得到较高的模拟精度. 有学者对随机格林法的地震动模拟进行一系列研究^[14-18]，在实测数据分析的基础上对Boore^[8]提出的随机格林函数法进行许多改进，提高格林函数法的模拟精度.

用随机格林函数法模拟地震动时，采用满足地震动傅里叶谱条件和持时包络曲线的随机波模拟子断层的地震动. 这种方法采用随机的相位谱，模拟结果很难保证地震动的方向性效应. 对此，有学者提出一些解决办法，如香川敬生^[14]提出选择位移波形相似的相位谱作为子断层的地震动，这种方法是从大量随机相位谱中挑选符合要求的相位谱，计算过程过于烦琐.

为了研究断层破裂方向和地震动传播特性对减隔震桥梁抗震反应的影响，本研究提出可以考虑地震动方向性效应的随机格林函数法相位谱设定方法，在此基础上利用北岭地震的断层参数，模拟不同震中距离的地震动作为地震输入，根据计算结果讨论场地与断层之间的相对位置对减隔震桥梁地震反应的影响，为桥梁减隔震设计的精细化分析提供一例参考.

随机格林函数法模拟地震动的基本原理是把断层分成若干个滑动量不等的子断层，各子断层按照断层破裂方向依次发生滑动，并向周围辐射地震能量，观测点的地震动U(t)为各子断层辐射而至的地震动u(t)的叠加结果. 许多学者在Boore^[8]的研究基础上进行一系列改进，已初步具备合理模拟设计地震动的能力. 本研究采用入倉孝次郎等^[15]根据相似比提出的地震动时程合成方法，即：

(1) $ U\left( t \right) = \sum\limits_{i = 1}^{N_{\rm{L}}} {\sum\limits_{j = 1}^{N_{\rm{W}}} {f\left( t \right)u\left( {t - {t_{ij}}} \right)} } , $

(2) $ \begin{split} f\left( t \right) = & \delta \left( t \right) + \frac{1}{{n'\left[ {1 - \exp \left( { - 1} \right)} \right]}} \times \\ & \sum\limits_{k = 1}^{\left( {N_{\rm{D}} - 1} \right)n'} {\left\{ {\exp \left[ { - \frac{{k - 1}}{{\left( {N_{\rm{D}} - 1} \right)n'}}} \right]\delta \left[ {t - \frac{{\left( {k - 1} \right)\tau }}{{\left( {N_{\rm{D}} - 1} \right)n'}}} \right]} \right\}}. \end{split} $

式中：t为时间，U(t)为大震的加速度时程，u(t)为小震的加速度时程，i为子断层沿断层长度方向的序号，j为子断层沿断层宽度方向的序号，k为细分后的小震沿断层滑动方向的序号，δ(t)为δ函数，τ为上升时间，N_L、N_W、N_D分别为断层在长、宽方向和断层滑动方向的分割数，n'为子断层滑动方向的细分数，t_ij为子断层ij地震动到场地的滞后时间

(3) $ {t_{ij}} = \frac{{{r_{ij}}}}{\beta } + \frac{{{\xi _{ij}}}}{{{V_{\rm{R}}}}} . $

式中：r_ij为子断层ij到场地的距离，ξ_ij为破裂开始点到子断层ij的距离，β为地震发生层的S波波速，V_R为破裂速度. 上升时间的计算公式为

(4) $ \tau = \alpha \frac{{{W_{ij}}}}{{{V_{\rm{R}}}}} . $

式中：α为0.25~0.60的经验系数，取值与断层类型有关；W_ij为对应子断层的宽度.

由于Boore^[8]提出的模拟方法不考虑地震动的近场和中场影响，影响长周期地震动模拟精度，为此本研究采用佐藤智美^[16]提出的考虑近场和中场影响的修正加速度傅里叶谱

(5) $ A\left( f \right) = {F_{\rm{s}}}S\left( f \right)P\left( f \right){N_{\theta \varphi }}\left( f \right) . $

式中：f为频率；θ、 $ \varphi $分别为震源球坐标系中的仰角、方位角；S(f)为震源特性；P(f)为传播特性； $N_{\theta \varphi }$为地震波传播近场和中场影响函数；F_s为自由表面放大系数，对于自由基岩表面，F_s=2.0.

(6) $ S\left( f \right) = C\frac{{{M_0}{{\left( {2{\text{π }}f} \right)}^2}}}{{1 + {{\left( {{f}/{{{f_{\rm{c}}}}}} \right)}^2}}}{\left[ {1 + {{\left( {{f}/{{{f_{\max }}}}} \right)}^n}} \right]^{ - {1}/{m}}} ， $

(7) $ C = \frac{{{R_{{{\theta }}\varphi }}P{\rm_{RTITN}}}}{{4{\text{π }}\rho {\beta ^3}}}{\left( {\frac{{\rho '\beta '}}{{\rho \beta }}} \right)^{0.5}} . $

式中：P_RTITN为地震能量在水平2方向的分散效果(0.71)；R_θφ为放射特性系数，取值与波的类型和频率有关，本研究取平均辐射系数(0.63)；ρ、ρ′分别为地震发生层和基岩层的密度；β′为基岩层的S波波速；M₀为子断层地震力矩；f_max为高频截断频率；m、n为常数，本研究取2；f_c为拐角频率，计算式为

(8) $ {f_{\rm{c}}} = 4.9 \times {10^6}\beta {\left( {{{\Delta \sigma }}/{{{M_0}}}} \right)^{{1}/{3}}} ， $

(9) $ \Delta \sigma = \frac{7}{{16}}\frac{{{M_0}}}{{{R^3}}} . $

式中：Δσ为应力降；R为等效圆形断层模型的半径.

(10) $ P\left( f \right) = \frac{1}{{{r_{ij}}}}\exp \; \left( {\frac{{ - {\text{π }}f{r_{ij}}}}{{Q\left( f \right)\beta }}} \right) . $

式中：Q(f)为与频率相关的传播路径衰减特性，具有地域性，本研究采用的计算式^[19]为

(11) $ Q\left( f \right) = 150{f^{0.5}}. $

近场和中场影响函数的计算式^[20]为

(12) $ \begin{split} {N_{\theta \varphi }} =& 6{\rm{i}}{\left( {\frac{\beta }{{{r_{ij}}\omega }}} \right)^3}\left\{ {1 - \exp \; \left[ { - {\rm{i}}\omega {r_{ij}}\left( {\frac{1}{{{\alpha _{\text{p}}}}} - \frac{1}{\beta }} \right)} \right]} \right\} - \\ & \begin{array}{*{20}{c}} \end{array} 6{\rm{i}}{\left( {\frac{\beta }{{{r_{ij}}\omega }}} \right)^2}\left\{ {1 - \frac{\beta }{{{\alpha _{\text{p}}}}}\exp \; \left[ { - {\rm{i}}\omega {r_{ij}}\left( {\frac{1}{{{\alpha _{\text{p}}}}} - \frac{1}{\beta }} \right)} \right]} \right\} + \\ & \begin{array}{*{20}{c}} \end{array} 2{\rm{i}}\left( {\frac{\beta }{{{r_{ij}}\omega }}} \right){\left( {\frac{\beta }{{{\alpha _{\text{p}}}}}} \right)^2}\exp \; \left[ { - {\rm{i}}\omega {r_{ij}}\left( {\frac{1}{{{\alpha _{\text{p}}}}} - \frac{1}{\beta }} \right)} \right] \hfill -\\ & 3{\rm{i}}\left( {\frac{\beta }{{{r_{ij}}\omega }}} \right) + 1 \hfill .\\ \end{split}$

式中：α_p为P波的波速，ω为圆频率.

随机格林函数法一般采用考虑时历特征的随机波模拟满足式(5)的地震动时程，计算过程如下：1）生成一组随机波；2）根据时程包络函数调整随机波在时间域的幅值；3）在频域调整各频率的幅值，让随机波的傅里叶谱与式(5)的结果一致. 这种方法存在地震动持时偏短、无法考虑地震动的方向性效应的问题.

本研究采用考虑地震动传播过程中弥散效应的经验时程包络函数^[17]，即

(13) $ {E_{^{_{{\rm{VL}}}}}}(t) = \left\{ \begin{gathered} {\left[ {{{(t - {t_{\rm{a}}})} \mathord{\left/ {\vphantom {{(t - {t_{\rm{a}}})} {({t_{\rm{b}}} - {t_{\rm{c}}})}}} \right. } {({t_{\rm{b}}} - {t_{\rm{c}}})}}} \right]^{\text{2}}},{\text{ }}{t_{\rm{a}}} \leqslant t < {t_{\rm{b}}}; \hfill \\ 1.0,{\text{ }}{t_{\rm{b}}} \leqslant t < {t_{\rm{c}}};{\text{ }} \hfill \\ \exp \; \left[ { - 10{{(t - {t_{\rm{c}}})} \mathord{\left/ {\vphantom {{(t - {t_{\rm{c}}})} {({t_{\rm{d}}} - {t_{\rm{c}}})}}} \right. } {({t_{\rm{d}}} - {t_{\rm{c}}})}}} \right]{\text{, }}{t_{\rm{c}}} \leqslant t \leqslant {t_{\rm{d}}}. \hfill \\ \end{gathered} \right. $

(14) $ \left. \begin{gathered} \lg ({t_{\rm{b}}} - {t_{\rm{a}}}) = 0.229M - 1.112, \hfill \\ \lg ({t_{\rm{c}}} - {t_{\rm{b}}}) = 0.433M - 1.936, \hfill \\ \lg ({t_{\rm{d}}} - {t_{\rm{c}}}) = 0.778\lg (r) - 0.340. \hfill \\ \end{gathered} \right\} $

式中：t_a为地震动开始时间，t_b为主震开始时间，t_c为主震结束时间，t_d为地震动结束时间，M为子断层的矩震级. 如图1所示为经验时程包络曲线形式.

本研究根据相位差分的分布模拟地震动时程^[18]. 子断层ij破裂引起的观测点地震加速度时程的傅里叶谱可以表示为

(15) $ {C_{ij,k}} = \frac{{{A_{ij,k}}}}{{{T_{ij,{{d}}}}}}(\cos\; {\varphi _{ij,k}} + {\rm{i}}\sin {\varphi _{ij,k}});\;k = 1,2, \cdot \cdot \cdot ,{N \mathord{\left/ {\vphantom {N 2}} \right. } 2} + 1 .$

式中：N为傅里叶变换的数据点数量；T_ij,d为子断层ij的地震动持时；φ_ij,k为子断层ij的相位；A_ij,k为子断层ij的傅里叶振幅，根据式(5)计算得到

(16) $ {A_{ij,k}} = {F_s}{S_{ij}}\left( {{f_k}} \right){P_{ij}}\left( {{f_k}} \right){N_{ij,\theta \varphi }}\left( {{f_k}} \right) . $

为了使式(15)的傅里叶系数经过逆变换后地震动时程与图1的持时曲线吻合，采用合理的相位谱. 由于持时包络曲线与相位差分谱相似，相位差分分布如图2(a)所示，如将相位差分谱等分成若干区段，可以得到如图2(b)所示的相位差分谱累计概率分布密度P_φ^[18]. 其中相位差分定义为

(17) $ \Delta {\phi _k} = {\phi _{k + 1}} - {\phi _k};\;k = 2,3 \cdot \cdot \cdot ,{N \mathord{\left/ {\vphantom {N 2}} \right. } 2}. $

根据傅里叶谱特性，φ₁=0, φ_N/2+1=0. 为了确定其他频率的相位，在0~1.0生成随机数λ_p，从图2(b)的累积概率分布密度中确定对应的相位差分Δφ_k，然后按照式(17)计算出相位φ_k+1. 当各频率的相位确定后，对式(15)的傅里叶系数进行逆变换，可以得到对应的地震动时程，结果满足式(16)的加速度傅里叶振幅要求，不需要迭代计算.

如果对每个子断层都按照上述方法生成地震动时程，由于子断层之间的相位没有规律，各子断层合成的地震动不能保证具有方向性. 实际上同一断层的相位具有相似性. 为此，本研究对第1个子断层外的其他子断层相位按式(18)计算，让小于拐角频率的各子断层相位具有相似性，且频率越小、相位越接近.

(18) $ {\phi _{ij,k}} = {\phi _{i{1},k}}(1 + {\mu _{ij,k}});\; k = 1,2, \cdot \cdot \cdot ,{N_{ij}}{,_{{\rm{fc}}}}. $

(19) $ {\mu _{ij,k}} = 0.1{({{{f_k}}}/{{{f_{ij,{\rm{c}}}}}})^2}(0.5 - \lambda );\; k = 1,2, \cdot \cdot \cdot ,{N_{ij}}{,_{{\rm{fc}}}}. $

式中：μ_ij,k为子断层相位相似性系数，f_ij,c为ij子断层的拐角频率，N_ij,fc为与拐角频率对应的傅里叶谱频域数据点. 当频率超过f_ij,c时，子断层之间的相位不考虑相似性.

子断层地震动模拟以及叠加过程的计算流程如图3所示.

以1994年美国北岭地震为例. 地震震级为M6.7，震中位于北纬34.211°，西经118.546°，震源深度为17.5 km. 根据Wald等^[21]的反演结果，断层走向为122°，倾角为40°，平均滑动角为101°. 断层破裂面沿走向的长度为18 km，下倾宽度为21 km，破裂面积为378 km². 断层顶部的深度为5.0 km，断层底部的深度为20.4 km. 应力降取5 MPa，地震矩估计为1.3×10¹⁹ N·m，断层面上的平均滑动量为101.85 cm，最大滑动量约为3 m. 震源附近区域的介质密度为2.8 g/cm³，震源附近区域的剪切波速为3.7 km/s. 断层面上的破裂传播速度约为震源附近区域剪切波速的0.8倍，取为3.0 km/s. 凹凸体数为2，面积分别为58.5、8.4 km²，滑动量分别为2.016、1.779 m；凹凸体地震矩为5.116×10¹⁸ N·m. 背景区域面积为310.8 km²，滑动量为0.661 8 m，地震矩为7.884×10¹⁸ N·m. 如图4所示为断层参数以及实测台站位置. 图中，OXYZ为整体坐标系，X为正北方向，Y为正东方向，Z为竖直向下方向，O点为断层面左上角点在地面的投影，oxy为断层面局部坐标系，o点为断层面左上角点，x为断层长度方向，y为断层宽度方向，Q点为断层面上破裂起始点的位置，ϕ_s为断层面走向，δ为断层面倾角，L为断层面长度，W为断层面宽度，d_t为断层顶部埋深.

如图5所示为北岭地震中3个台站模拟的地震动加速度时程和反应谱与相应实际记录的对比，计算时相位差分分32等分计算累计概率密度分布. 图中，t为时间，a为加速度，a_max为最大加速度，T为周期，S_a为加速度反应谱，H₁、H₂分别为水平相互垂直的2个方向. 结果表明，无论是峰值加速度还是宽频域的反应谱，均有比较好的精度. 由于震源参数、传播路径以及场地条件的不确定性，以及计算方法的近似性，计算结果与实测结果多少存在差异. 实际地震记录取自文献[22].

为了验证模拟地震动的速度脉冲参数准确性，表1给出位于破裂方向上台站LV3、LV1地震动的PGV/PGA. 可以发现，本研究模拟的地震动PGV/PGA与实测记录比较接近.

以如图6所示的铁路减隔震桥梁^[23]为例. 桥梁为9 m×32 m的预应力混凝土简支箱结构，下部结构采用U型桥台双线圆端型空心桥墩. 桥墩编号从左至右依次为P1~P8，梁端桥台分别为A1、A2. 桥墩高度依次为24、35.5、41、42、33、31.5、22、10.5 m. 采用E型钢阻尼作为减震元件. 支座编号从左至右依次为B1~B18. 如图7所示为桥梁−轨道一体化计算模型，取道床阻力系数为15 kN/m.

断层条件仍采用北岭地震数据. 如图8所示选取4个方位的观测点的地震动作为顺桥向输入条件，分析方向性效应、上盘效应对减隔震桥梁地震反应的影响. 图8(a)为设定的4个观测点与断层的位置关系. 震中距离设为15、30 km共2种. 其中观测点1、2分别位于地震动的破裂方向以及反方向，3、4位于断层的2个侧面，其中测点4位于上盘一侧. 为了方便表述，将2种震中距离的观测点分别编号为15-1、15-2、15-3、15-4、30-1、30-2、30-3、30-4。图8(b)、(c)为模拟地震动时程以及反应谱. 不难看出，模拟的地震动能较好地反映出地震动的方向性和上盘效应，观测点1的地震动明显大于其他观测点，位于上盘的观测点4大于下盘的观测点3. 当震中距分别为15、30 km时，PGV/PGA分别为0.13、0.11，均小于0.2，尚未达到近场地震动反映速度脉冲效应的程度.

如图9所示为震中距离分别为15 、30 km时，桥墩P4的墩顶位移时程以及P5桥墩E型钢阻尼器的荷载−位移滞回曲线. 图中，d为墩顶位移，F为支座反力，d_max为墩顶最大位移. 由图可知，当桥梁位于具有方向性效应的观测点1时，墩顶位移以及阻尼器的地震反应明显大于其他观测点，位于上盘的观测点大于等震中距离的下盘观测点. 表明仅根据震中距离不考虑震源条件的结构抗震设计方法存在一定的缺陷，存在预测的结构地震反应过大或者偏小的问题. 还可以发现，在震源附近虽然方向性效应导致不同观测点地震动的峰值有明显的差异，但是这些观测点的地震反应谱特性相似，且桥梁的地震反应规律相同. 如震中距30 km破裂点方向上的观测点峰值加速度与震中距15 km断层侧面点的观测点峰值加速度相近，绝对值分别为3.14、3.02 m/s²，2种波输入时计算得到桥梁P4墩的墩顶最大位移绝对值分别为25.6、26.2 mm，P5墩E型钢最大位移绝对值为42.3、45.3 mm. 即峰值加速度相同的地震动，其地震反应比较接近，没有出现由于桥位位置不同导致反应明显差异的情况. 这是因为这些台站的地震动速度脉冲指标尚未达到影响减隔震桥梁地震反应的程度，加速度特性仍然是控制结构地震反应的主要因素.

为了进一步比较不同地震动加速度时程对减隔震桥梁地震反应的影响，选用常用的地震动时程El Centro波作为输入进行分析，峰值加速度分别调整到0.15g、0.30g 、0.6g，将地震动时程进行适当调整，使其反应谱调整至与观测点1地震动的反应谱相同. 如图10为调整后的El Centro地震动加速度时程.

表2为桥墩P4墩顶的最大地震位移绝对值 $|{{d}}^{{\rm{P4}}}_{ {\rm{max}} }|$以及P5的E型钢阻尼器塑性率η^P5计算结果. 结果表明，用相同的加速度反应谱输入时，El Centro地震动输入时的结构地震反应与模拟地震动时程的结果有一定的差异，当PGA=0.60g时，模拟地震动输入时的结构地震反应明显大于El Centro地震动，但当PGA=0.30g时，模拟地震动输入时的结构地震反应却比El Centro小. 因此，减隔震桥梁的地震反应受到地震动输入影响大，有必要在设计时输入合理的地震动加速度时程。考虑断层破裂过程模拟的地震动，反映断层参数的影响.

（1）用相位差分谱模拟随机格林函数法的地震动加速度时程，计算过程方便，模拟的地震动时程能够反映持时包络曲线的基本特征.

（2）采用考虑近场项、中间项的随机格林函数法模拟的地震动，在较宽的频域内反应谱与实测结果的基本接近，且PGV/PGA接近.

（3）考虑子断层之间的相位谱在低频域的相似性，随机格林函数法可以反映断层的方向性效应。

（4）在近断层地震中，当PGV/PGA较小时，反应谱特性基本相似，减隔震桥梁的地震反应主要与地震动的强度有关，与断层−桥位的方位基本无关.

（5）减隔震桥梁的地震反应受到输入地震动加速度时程影响大，即使相同弹性反应谱的地震动输入，也会得到不相同的地震反应.

（6）本研究采用相位谱模拟地震动加速度时程，对考虑地震破裂过程的减隔震桥梁地震反应进行了探索，关于考虑子断层相位相似性的公式还应从理论上进一步完善和验证.