0 引言

大震的破裂尺度一般为几十千米到几百千米, 其震源不是一个点而是一个面。如四川汶川8.0级地震破裂长度约200km, 其烈度等震线呈狭长椭圆状, 具体表现为沿地震破裂方向的烈度衰减缓慢, 等震线的长轴尺度远远超过了点源烈度衰减关系所预测的尺度；垂直于地震破裂方向的烈度衰减较快, 等震线的短轴尺度与烈度衰减关系所预测的尺度基本一致(帅向华等, 2020)。因此, 基于线源的烈度衰减模型能够更准确地评估烈度分布, 对地震破裂方向和破裂方式的快速准确判定十分重要。在目前常用的地震烈度评估模型中，地震破裂方向的判定方法相对简单, 即将震中位置与地震活动断裂带数据进行叠加分析, 计算震中到周围断层的最小距离, 并选择距离最近的断层走向作为地震破裂方向(地震影响场长轴方向), 且破裂方式均为对称双侧破裂。这种方法的优势在于时效性高, 可快速获得评估结果。但该方法在实际工作中存在一定的局限性：一是地震活动断裂带数据不完备, 尽管经过地质专家多年的调查, 已经积累了大量的断裂带数据, 但仍然不排除有未发现的隐伏断裂, 如2014年8月3日云南鲁甸6.5级地震的发震断层实际为NW向的包谷垴-小河断裂, 而当时调查记录的距离震中最近的为NE向的昭通-鲁甸断裂；二是地震有时会发生在2条或多条断裂交汇的地方, 此时震中到几条断裂的距离十分接近, 距离最近的断裂不一定为发震断层；三是大震的破裂方向往往是单侧破裂或者不对称的双侧破裂, 如四川汶川8.0级地震以NE向为主、四川芦山7.0级地震以SW向为主。因此, 迫切需要一种快速且可靠性较高的地震破裂方向的判定方法, 提高烈度速报的效率和精度, 为地震应急救援及灾害损失快速评估提供更加快速、准确的决策支撑和服务保障。

地震断层破裂传播的方向对地震波(体波和面波)的影响可以用有限移动源模式下的地震波谱理论来解释。有限移动源模式是指地震断层面上的各点不可能同时发生错动, 通常为先从某一点(或部分)开始破裂, 然后以有限速度向断层的其他部分传播。因而在地震记录中必定带有某些震源的相关信息, 如震源破裂速度等其他参数的信息。相关研究结果表明, 震源的有限性和有限的破裂速度对地震波的影响主要体现在2个方面：一方面, 地震波的振幅在方位上受到有效的调制, 也就是所谓的“方位调制”作用, 即地震波振幅在沿着破裂的传播方向上加强, 而在相反方向上减弱；另一方面, 地震波卓越周期随方位有类似于多普勒效应的变化(胡进军等, 2011)。因此, 根据体波振幅和卓越周期在不同方位上的变化即可反推出地震断层的破裂方向。

趋势面分析是一种多元统计分析方法, 能有效地分析某一种属性数据在空间上的分布规律与变化趋势。该方法最初应用于地理科学中, 利用数学趋势面模拟地理系统要素在空间上的分布规律及变化趋势(王江萍等, 2009), 后来其被广泛地应用于模拟环境、资源、疾病、人口、农业、经济等要素在空间上的分布规律(Hota, 2014)。

本文根据有限移动源理论, 提出一种利用趋势面分析和震中一定范围内台站初始P波地震动参数判断地震断层破裂方向的方法, 并基于断层距模型和最小二乘法统计回归确定地震断层破裂方式。通过实际震例, 对比分析根据本文所提方法评估的地震断层破裂方向与实际地震影响场长轴方向的吻合度, 评估采用该方法进行地震烈度速报的可靠性, 为利用地震监测台站实时数据进行地震断层破裂方向判定提供一种新思路, 为地震烈度速报结果的动态修正提供参考。

1 地震断层破裂方向研究现状

确定地震断层破裂的方向主要有现场地质探查、大地测量观测以及基于地震学资料的反演分析3大类方法(何骁慧等, 2015)。现场地质勘察方法为地质学专家在震区考察识别地震破裂带、测定断层位置及走向, 是判定断层面直接有效的方法, 其得到的结果最精细, 但是所需时间较长, 不适合地震烈度快速评估需求。大地测量观测方法采用合成孔径雷达干涉测量技术(InSAR)测定地表形变, 进而测定断层破裂方向, 该方法仅在地表裸露、植被不茂密的地区较为有效, 且因卫星重复轨道数据往往无法快速获取, 不能用来快速(几个小时或者当天)测定破裂断层。地震学资料反演的方法主要是根据余震的空间分布、强震动台观测的地震动信息(峰值速度PGV、峰值加速度PGA等)或多个测震台上记录到的波形信息推断地震断层破裂的方向, 从而确定发震断层。其中余震分布法一般需要震后2h内的余震精定位信息, 时效性较差, 而我国强震动台站数量较少的现状会影响破裂方向评估的准确度, 因此, 基于测震地震波形判定地震断层破裂方向亟待研究。

经验格林函数法是研究地震破裂方向性的常用方法之一, 该方法使用小震作为经验格林函数校正地震波在三维地球介质中的传播效应, 从而对震源的特性进行更准确的估计(Mozaffari et al，1998)。Tan等(2010)在此基础上, 利用不同方位台站记录P波的形状和振幅信息, 发展了正演模拟技术, 从而测定破裂方向。对于大震(M>7), 简单的破裂方向参数已不足以描述破裂的复杂性, 一般使用有限断层反演方法研究其破裂过程的细节(王卫民等, 2008、2013)。另外, 通过地震质心位置与破裂起始位置之间的差异, 也可测定主震的破裂方向以及发震断层的破裂方向性(秦刘冰等, 2014)。上述研究中基于震源参数反演方法所需处理的资料在主震发生后很难快速得到。冯蔚等(2015)利用测震台站波形数据, 计算LQT坐标系下得到的S波振幅值与单力偶源S波辐射图案匹配, 判断鲁甸地震的发震断层方向与调查结果吻合, 但是对于大震, 近场测震台站S波记录受限幅影响较大(冯蔚等, 2015)。在地震破裂传播方向上, 地震波的周期变短；在其相反方向上, 地震波的周期变长。对一些较大地震, 在宽频数字化记录中能够比较清楚地观测到地震多普勒效应(Douglas et al，1981), 也有学者研究了利用P波的多普勒效应来判断地震破裂的方向(周云好等, 2002)。

一次大震的断层破裂过程非常复杂, 会发生多次破裂, 主震的首次破裂一般释放的能量最大, 对震源参数的确定起主导作用。随着地震预警技术研究和发展, 许多学者发现初始P波能够较好地反映主震的震源特征。研究表明台站记录的P波前3s最大位移与相应台站的PGV具有较好的相关性(Wu et al，2005)。Colombelli等(2015)在此基础上利用日本KiK-net台网记录的本州岛海域9.0级地震及其余震波形记录的P波初始3s幅值来估算台站的PGV, 再由PGV与修正的麦卡利烈度IMM的关系来计算仪器烈度, 最后通过内插来评估地震烈度分布。

2 数据准备

本文选取了2008—2014年8个破坏性地震作为实验震例(表 1), 其中M_S8.0地震1个, M_S7.0地震1个, M_S6.0~7.0地震6个, 震源深度从最浅5km到最深的20km, 分布在四川、云南、甘肃和新疆4个多震省份。

实验数据选取了以上震例中距震中300km范围内的测震台站和强震动台站连续波形记录, 并对波形记录进行预处理(图 1), 经过基线校正和去线性化处理, 对波形记录的误差进行修正, 采用了0.075~10Hz的中通Butterworth滤波器进行滤波, 然后利用STA/LTA方法进行P波震相的自动拾取, 以1s为间隔截取2~10s的波形记录, 分别对记录进行积分、水平向合成、三分量合成，即

$ v({t_i}) = \sqrt {v{{({t_i})}_{{\rm{E}} - {\rm{W}}}}^2 + v{{({t_i})}_{{\rm{N}} - {\rm{S}}}}^2 + v{{({t_i})}_{{\rm{U}} - {\rm{D}}}}^2} $

(1)

其中, v(t_i)为三分量合成速度, v(t_i)_E-W为EW向速度, v(t_i)_N-S为SN向速度, v(t_i)_U-D为垂直向速度。

之后取其极值得到不同时段垂直向、水平向和三分量合成的峰值速度(Peak Velocity, 以下简称PV)、峰值位移(Peak Displacement, 以下简称PD)以及对不同时段的记录进行FFT变换得到卓越周期(Dominant Period, 以下简称DP)。最后将计算数据写入表格, 形成数据文件, 共得到240组台站记录。

3 方法 3.1 破裂方向分析方法

破裂方向分析采用趋势面分析法。趋势面分析法是研究一定区域内空间数据发展或演变趋势的一种方法, 该方法的实质是数据拟合(查文婷等, 2015)。将局部区域内某一事件的观测值分为趋势值和残差值, 趋势值反应了该事件空间要素的宏观分布规律, 可以排除一些偶然因素引起的变异影响。自变量为观测点的经纬度坐标或平面坐标, 因变量为观测值, 用最小二乘法获得残差平方和最小估计值, 以反映观测值变化的整体趋势。本文采用的是多项式拟合, 即利用线性模型对某一变量的观察值进行多元回归拟合, 以产生该区域的趋势面。趋势面分析的多项式函数表达为

$ Z = {{\rm{A}}_1} + {{\rm{A}}_2}x + {{\rm{A}}_3}y $

(2)

其中, x、y为观测点的平面横、纵坐标, Z为观测值。

3.2 破裂方式分析方法

对于本文研究的6级以上强烈地震, 震源可以近似为一条地震沿断层破裂的线段。破裂断层一般分布于震中的一侧或两侧, 取决于地震破裂方式。如图 2所示, 设破裂断层一端到震中的长度ao占断层总长度ab的比例为per。在震中位置、破裂方向和破裂长度确定的情况下, 变换per比例即可得到所有可能破裂方式的断层分布。

观测点的震中距应取其到破裂断层线段最近的距离D(图 3)。当观测点位于破裂断层两端之外时(如图 3中a点、d点), 取其到断层两个端点的距离；当观测点位于断层两个端点之间时(如图 3中b点、c点), 取其到断层的垂直距离。

本文采用常用的地震动衰减关系模型为(Kanai, 1961)

$ {\rm{lg}}Y = {{\rm{C}}_1} + {{\rm{C}}_2}M + {{\rm{C}}_3}{\rm{lg}}(R + {R_0}) $

(3)

其中, Y为地震动参数, 本文选取的是三分量合成的峰值速度(PGV)；M为震级；R为震中距, 本文采用断层距表示；R₀控制近场Y为一个有限值, 由于近场地震动的饱和特性, 可以假定震中一定范围内的地震动饱和, 所以本文设定R₀为10km, 即距震中10km内地震动饱和。

将观测点的断层距与PGV代入式(3), 利用最小二乘法进行回归分析, 每种断层破裂方式可以拟合出一个结果, 这些结果对应的标准差越接近0, 说明拟合的误差越小, 也就是说最小标准差对应的破裂方式与观测值的匹配度最高。本文设定per比例为0到100%之间, 变换步长为5%, 通过迭代回归分析21种破裂方式的衰减模型, 以拟合标准差最小的per比例作为断层破裂方式评估结果。

4 结果及分析 4.1 破裂方向分析结果

本文利用ArcGIS软件对一定分布范围内台站的初始P波参数进行插值分析以获得地震断层破裂的趋势方向, 评估流程如图 4所示。其中, 在选取插值范围时, 应保证震中的各个方向上均选取一定数量台站, 且台站分布较均匀。

破裂方向分析以四川芦山7.0级、云南景谷6.6级和四川汶川8.0级地震为例, 采用了距震中200km范围内的测震台站波形记录, 选取了P波初始3s的最大位移和初始5s的卓越周期, 再利用反向距离加权和趋势面插值方法分析地震影响场长轴方向, 研究结果表明：芦山地震P波幅值和卓越周期的反向距离插值IDW结果呈现出一定的破裂方向性, P波初始3s峰值沿发震断层NE方向较低, SW方向较高, 5s内的卓越周期沿发震断层NE方向较高, SW方向较低, 表现出芦山地震不对称破裂方式的特点, 2种参数的趋势面插值结果得到的趋势方向均与实际调查烈度影响场的长轴方向基本一致(图 5、图 6)；景谷地震P波幅值的插值结果与该地震较对称的破裂方式特点一致, 即P波初始3s峰值沿发震断层的2个方向上没有明显差别, 但是5s内的卓越周期的插值结果呈现出延发震断层SE向的破裂趋势, 2种参数的趋势面插值结果得到的趋势方向均与实际调查烈度影响场的长轴方向基本一致(图 7、图 8)；汶川地震P波幅值的IDW插值结果具有比较明显的破裂方向性, 即P波初始3s幅值沿发震断层NE方向较高, SW方向较低, 而且有显著的不对称破裂特点, 趋势面插值得到的趋势方向与实际调查烈度影响场的长轴方向一致(图 9), 汶川地震5s内的卓越周期不具有显著方向性。

在实验过程中发现, 初始时刻和P波参数的选择对分析结果有较大影响。因此, 本文利用实验震例数据, 分析不同初始时段PD、PV和DP进行趋势面分析的可靠性。具体方法如下, 分别利用初始2~10s的PD、PV和DP进行趋势面分析, 将分析得到的趋势方向与影响场方向进行比较, 完全一致计3分, 基本一致计2分, 不一致计0分, 以打分之和代表某一时段对应参数预测的准确率。当趋势方向与影响场方向一致时可判断为正确, 正确震例占所有震例的比例即为正确率。

分析结果表明，PV的准确率总体最高, 其次是PD, DP最低, 即参数的可靠性PV>PD>DP (图 10)。PD在4s和5s时正确率最高(83%)；PV在4s、9s和10s时正确率最高, 均达到了100%；DP在7s时正确率最高, 也达到了100%。所有参数的总体正确率达到70%以上；DP在6s前的正确率较低, 原因可能是近场台站密度较低, 6s内的卓越周期在近场差异性不够显著。

4.2 破裂方式分析结果

破裂方式分析以汶川8.0级地震为例, 利用震中300km范围内台站初始P波3s位移峰值的趋势面分析得出断层破裂走向为NE∠50°, 破裂长度计算采用Wells等(1994)关于破裂长度L与震级M之间的统计关系

$ \lg L = 0.69M - 3.22 $

(4)

汶川地震断层破裂长度由式(4)计算约为200km。在确定方向和破裂长度情况下, 变换per比例可得到21个可能的断层分布。根据线源断层距衰减模型(式(3)), 对观测台站lgPGV拟合结果的标准差如图 11所示, 最小标准差对应的per为0, 即评估破裂方式为NE∠50°的单侧破裂。

将本文利用台站观测数据和模型评估计算得到的汶川地震破裂断层分布绘制到地图上, 与汶川地震现场调查烈度分布进行对比(图 12)。评估破裂断层从震中汶川县沿NE向经北川县延伸至青川县, 评估的断层破裂方向和破裂方式与现场实际调查结果基本一致。

5 结论

在地震烈度速报中, 对地震断层破裂方向效应的考虑十分重要, 如何快速评定地震断层破裂方向和破裂方式成为相关研究的关注点。本文提出了一种利用趋势面分析法和台站观测数据快速评估地震破裂方向和破裂方式的方法, 并基于近年来我国西部发生的一些强震实例分析了该方法的可靠性和实用性。主要研究结果如下：

(1) 提出的基于趋势面分析和观测台站初始P波的地震断层破裂方向快速判断方法具有较高的可行性, 对于不同区域、震源深度及破裂类型的中强地震均具有较好的适应性。但方法的准确率受地震动参数(速度峰值、位移峰值和卓越周期)的选取、P波初始时段选取及台站选择等影响, 对于震中附近台站非常稀疏的地震无法适用。

(2) 在初始P波信息处理中, PD的最优初始时段选择为3~5s和9~10s, PV的最优初始时段选择为4~5s和9~10s, DP的最优初始时段选择为6s以后；该方法使用的地震动参数可靠性为PV>PD>DP, 但总体正确率较高, 达到了70%以上。

(3) 提出利用台站观测数据和最小二乘法拟合不同破裂方式的PGV衰减关系, 以拟合结果最优的破裂方式来确定地震破裂断层分布的方法具有可行性, 在震后能够利用实时获取的台站数据快速评估破裂断层分布, 能够较大地提高地震破裂断层评估的效率。

(4) 实验震例中7级以上地震较少, 因此该方法对于大震破裂方向判定的可靠性还有待检验。另外, 本文选取的测震台站分布密度较低, 需要进一步利用密度较高的地震监测台网来进行记录, 以检验该方法的适用性。

(5) 本文快速判定地震破裂方向的整个处理过程包括波形数据预处理、P波震相拾取、一定时窗内P波峰值速度或峰值位移参数提取, 然后通过趋势面分析法评估地震破裂方向, 再利用最小二乘法分析地震破裂方式, 最终得到破裂断层分布。整个过程可以通过编程自动计算完成, 能够构成一个系统，处理计算时间不超过1min, 震后10min内可以产出判定结果。