引言

地震预警系统具有减轻地震灾害的能力，在过去的数十年里，地震预警已经逐渐发展成为减轻地震灾害的一个重要手段（Nakamura，1988；Allen，Kanamori，2003；Zollo et al，2010 ；金星等，2012）。世界上许多地区建设了由实时监测地震的台站和中心处理系统构成的地震预警系统，包括日本（Odaka，2003；Horiuchi et al，2005 ；Kodera et al，2016 ）、美国加州（Allen，Kanamori，2003；Allen et al，2009 ；Böse et al，2009a ）、吉尔吉斯斯坦（Pittore et al，2014 ）、墨西哥（Espinosa-Aranda et al，1995 ，2009）、意大利（Zollo et al，2006 ）、罗马尼亚（Böse et al，2009b ）和伊比利亚半岛西南部（Pazos et al，2015 ）。我国北京、福建和台湾地区也陆续建成了地震预警系统，开展了一系列测试工作（Wu et al，1999 ，2013；马强，2008；Peng et al，2011 ；Chen，2015）。目前我国的国家地震烈度速报与预警工程正处于实施阶段，预计2020年投入测试工作，有望在全国范围内提供地震预警与烈度速报服务。

大地震往往沿发震断层主破裂方向产生较长距离的破裂，如M_W7.0约60 km，M_W7.5约130 km，M_W8.0约300 km （Wells，Coppersmith，1994）。地震预警系统在对地震的震中、震级等参数进行估测时，通常将地震震源视为一个点源，这样的点源假设在实际应用中可以较好地表示较近距离范围内中小震的震源，但对于罕遇的大震（M＞7.5），则会出现明显的偏差。2011年日本东北（Tohoku） M_W9.0地震中，在仅使用P波信息的条件下，地震预警系统给出的估测震级仅为M5.7，东京地区的烈度估测结果为日本气象厅地震烈度（Japan Meteorological Agency，简写为JMA）4度，由于预警系统未考虑断层破裂的长度，包括东京在内的日本关东部分地区均出现了烈度被低估的现象，东京地区的实际烈度达到了JMA5度（Hoshiba et al，2011 ）。地震预警系统在面对这样的大震时，如果能够迅速地估测出地震的破裂长度，则可通过估测的破裂长度修正震级和地震动参数，从而降低低估的可能性。

大地震破裂长度的快速估测近年来一直是地震预警系统研究的重要方向，而且相关研究也取得了一些突破。土耳其伊斯坦布尔的预警系统建设中采用PreSEIS方法（Böse，2006），该方法基于神经网络方法，将点源方法扩展到随机有限断层模型中，通过对历史有限断层模型数据集的训练，学习得到断层位置、破裂长度、地震震级等信息。这种方法能够在地震发生时给出地震破裂长度的估测值，且估测速度与现有的原地预警方法接近。Böse等（2012）开发了有限断层破裂探测器“FinDer”算法，该方法基于密集的地震台网和一组预先计算好的断层模板，在地震发生时实时测定地面运动的振幅，采用图像识别技术自动实时计算断层破裂（假设为线源）在地表投影的中心位置、长度和走向；Zhang等（2014）研究了利用近场强震动数据进行地震破裂过程实时成像的方法。近年来，随着近场台阵反演震源破裂过程的研究取得较好成果（Fletcher et al，2006 ；Vallée et al，2008 ；Meng et al，2012 ），研究人员开始尝试利用近场台阵来为地震预警系统提供信息。Meng等（2014）提出了利用近断层密集台阵波束形成技术反演震源方位角的方法来实时估测震源的破裂长度，用美国的UPSAR强震动台阵记录对多次地震进行破裂长度的估测，取得了较好的结果。该方法需要小孔径地震台阵提供精确的走时，UPSAR台阵能够较好地满足这一要求，但我国现有的一些强震动观测台阵（例如，在汶川地震中的四川自贡强震台阵）中各观测子台的观测仪器均为独立授时，各子台记录可能存在时钟同步误差问题，直接应用该算法还需要寻找适合每个台站的延时确定方法。因此，本文拟基于波束形成技术开发一套通过反演破裂主能量源方位角的震源破裂长度快速估算算法，同时提出调整各台站相位的方法，即在每个台站的延时无法确定的情况下，能够近似地得到每个台站之间的相位差；并在此基础上对2008年汶川地震四川自贡强震台阵记录到的强地震动记录进行分析，以探讨该方法对汶川地震的断层破裂长度进行实时估测的可行性。

1.   方法

利用近断层密集台阵记录反演破裂主能量源方位角的方法来实时估测震源的破裂长度，涉及波束形成技术，而其中的密集台阵中台站之间时间延迟（也就是相位差）的确定及其精度则是关键所在。

1.1   波束形成技术

波束形成技术是阵列信号处理的重要部分，该方法通过将各阵元输出进行加权求和，以达到增强期望信号、抑制干扰的目的。利用阵列信号处理技术分析相干信号的相位，可以确定波在阵列上的传播方向（张小飞等，2015），因此，可以通过密集地震台阵信号确定某一时刻地震波的传播方向，进而确定产生这部分地震波的破裂源方向。破裂源方向射线与预定义的线性断层的交点则可视为破裂的主能量源位置，也可将其视为地震的瞬时震中。基于这些瞬时震中即可估测断层破裂的尺度，即将每一时刻计算出的瞬时震中之间的最大距离视为断层破裂长度。

在波束形成过程中采用每两组波形之间的互相关系数进行叠加来代替原始波形的叠加，能够减少地震波在地壳浅部环境传播过程中的散射和多路径所造成的不利影响（Borcea et al，2002 ；Fletcher et al，2006 ）。因此，对于一个选定的时窗，将每两组地震动记录之间的互相关系数c_ij（θ）定义为

式中，x_i为第i个台站记录到的地震信号，t为时间，τ_ij（θ）为台站i与台站j之间的时间延迟，

式中：θ为慢速矢量，对于单频波，其慢速矢量定义为大小为波速的倒数、方向与单频波传播方向相同的矢量；δt_i和δt_j分别为台站i和j本身的记录延迟。将求得的τ_ij（θ）舍入到最近的采样点上用于计算c_ij（θ）。

这样，按照式（1），在每一个时间窗下，求出每两组地震动记录之间的互相关系数c_ij（θ）并求和得到互相关系数的和c_ij ^s （θ）。由于慢速矢量包含着地震波传播视速度和反向方位角的信息，c_ij ^s （θ）是慢速矢量θ的函数，因此对于每一个时窗而言，当c_ij ^s （θ）的极值大于某个规定的阈值时，就认为此时θ所对应的方位角是地震主要辐射体所在的方向。Fletcher等（2006），Vallée等（2008）和Meng等（2012）均利用这种方法来推断破裂长度、断层分支及破裂速度，但其中涉及较多的人为推断，在实时应用中难以实现。

1.2   改进的相位调整方法

波束形成方法需要台阵中每个台站准确的延时数据δt_i。然而，现有的强震动观测台阵的强震动记录可能存在丢头或时钟同步误差较大等问题，因此在对各地震动时程进行分析之前，需要确定台阵中每个台站相对于标准台站的相位差。基于这个需求，Fletcher等（2006）曾提出利用大量余震、小震数据统计分析计算每个台站延时的方法，并取得了良好的效果；但对于本文研究所用的已获取汶川地震主震记录的自贡强震动观测台阵而言，能够获得的余震、小震记录少，无法通过该方法来计算台阵的延时。另外，Meng等（2012）曾利用10 s时窗下的相关性分析方法进行相位调整，但该方法旨在分析断层的破裂过程，10 s时间窗的设置对于分秒必争的地震预警系统显得尤为过长。

为解决上述台阵各台站延时确定的问题，本文提出相位调整的方法，具体步骤如下：

1） 利用Allen （1978，1982）的P波捡拾方法，对地震动时程的垂向分量进行首波初步捡拾；

2） 采用Akaike （1973）提出的赤池信息准则（Akaike information criterion，简写为AIC）进行首波的精确捡拾；

3） 给每一条地震动时程添加一个3 s的时窗，时窗起始点为首波的捡拾点；

4） 进行0.2—0.7 Hz带通滤波；

5） 利用相关性分析方法剔除丢头的数据；

6） 调整每个台站间水平向、垂向空间位置所造成的相位差，其中：水平向空间位置造成的相位差近似为台站空间位置矢量在点源震中到台阵的方位角上的投影长度除以莫霍面上的纵波波速；垂向空间位置造成的相位差近似为台站的高程差除以上地壳的纵波波速。

经过上述步骤，能够剔除丢头的地震动时程，获取满足波束形成技术的地震动时程。

图1给出了本文方法的技术路线，可以看出，该算法主要分为两个部分：第一部分为数据筛选与相位调整，该部分依照相位调整方法估算各台站相对于参考台站的相位差，并根据相位差调整每组数据的相位；第二部分为波束形成方法，按照1.1节所介绍的方法对断层破裂长度进行估算。

图  1  本文的技术路线方案

Figure  1.  The technical method of this paper

下载: 全尺寸图片 幻灯片

2.   基于汶川地震强震动台阵记录的应用分析

自贡地形影响强震动观测台阵是“十五”强震动台网系统建设项目的组成部分，于2007年建成并实施观测运行。该台阵位于我国四川省自贡市西山公园内的山体上，由8个台站组成，每个台站均配置了ETNA型记录器和ES-T型加速度计的数字强震动仪（杨宇等，2011；唐晖等，2012）。山顶台站（6号）与山底台站（0号或1号）的最大高程差达到72 m，4号台站与5号台站之间的高程差最小，为2 m。0号台站位于山底土层之上，1号台站位于山底基岩上，其它6个台站位于山脊不同高度的侏罗纪基岩上，各台站之间的空间相对距离处于52—400 m之间（阿布都瓦里斯，2013）。

图  2  自贡强震台阵记录到的汶川地震三分量加速度时程

Figure  2.  The three-component acceleration time histories of Wenchuan earthquake recorded by Zigong strong motion array

下载: 全尺寸图片 幻灯片

自贡台阵位于汶川地震震中东南约200 km，在汶川地震中，该台阵获得了良好的主震加速度时程，各台站的加速度时程如图2所示，可见记录到的地震动时程的近震特征明显。

利用地震台阵估测地震主能量源反向方位角的方法是基于地震台阵信号的相位和相干性来追踪产生最强辐射源的区域。本文将2号台站作为参考台站来讨论台站之间的相位差。从图2中可以看出，1号台站的地震动时程相比其它台站数据存在明显的相位前移。为了筛选可用的地震动记录并对齐首波，需要估算地震波的到时。采用吴建平等（2009）研究中的汶川地区地下速度结构（图3b）以及刘瑞丰等（2014）的震相标识方法，粗略地估计出地震动记录各标志震相的到时；然后根据到时估算出P波与S波的到时差约为27 s。图3a为1号台站地震动时程的S波标志震相标记，可见Sn，Sb和Sg震相大概出现在15—21 s，考虑到估算的P波与S波的到时差，推断1号台站的地震动时程出现了丢头的现象，再加上波束形成方法需要对齐初至的Pn波相位，所以1号台站的地震动时程不能满足本研究需求。

台阵中的每个台站均存在各自的延时，其延时δt可以由大量的余震记录通过统计方法获得，由此精确地校正地震动时程的初始相位。由于可获得的自贡台阵数据量较少，我们无法通过余震记录的统计方法获得每个台站的δt，所以首先对齐首波，再根据各台站的空间相对位置确定每个台站与参考台站之间的相位差。本文采用P波捡拾方法（Allen，1978，1982），对地震动时程的垂向分量进行首波的初步捡拾，然后采用AIC方法（Akaike，1973）进行首波的精确捡拾，以此来对齐时程。

从各个台站之间Pn波到达后的初始10 s时间窗下波形的相关性（图3c）来看：由于1号台站记录到的地震动时程存在丢头的现象，其与其它台站之间的波形相关性非常小；而土层上的0号台站由于受到放大效应的影响，其地震动时程与除1号外其它6组地震动时程的相关性较差。在台阵中台站数量足够大的情况下，阵列数据处理中一般选择相关系数较高的地震动时程（Meng et al，2012 ）。在实时应用中，本文通过设置相关性系数的阈值，实现了自动择取相关性较好的、适合估测断层长度的数据的算法。考虑到本文可用的台站数据数量较少，计算中将阈值设置为0.3，剔除1号台站数据，选用了其余7个台站的数据。

图  3  强震动记录的筛选过程示意图

（a） 1号台站前60 s三分量强震动记录及其S波标志震相；（b） 本文采用的汶川地区地下速度结构；（c） 自贡台阵记录到的汶川地震各台站加速度时程之间的相关性

Figure  3.  Diagram of the screening process for strong ground motion records

（a） The pre 60 s three-component strong motion records of the station No.1 and the seismic phase of the S wave；（b） Underground velocity structure of seismic wave in Wenchuan area；（c） Correlation between acceleration time histories of each station recorded by Zigong array during Wenchuan earthquake

下载: 全尺寸图片 幻灯片

由于台站的空间位置不同，对齐首波的时程并不能反映各台站之间的真实相位差，因此本文依据震中到台阵的方位角和各台站的空间位置计算各台站与参考台站之间的相位差，并根据相位差对各台站的波形相位进行调整。水平向的相位差由各台站的经纬度坐标、点源震中到台阵的方位角以及莫霍面的波速决定；垂向的相位差理论上可以通过多次余震的数据统计得到（Fletcher et al，2006 ），但由于自贡台阵缺少相关数据，故本文选用台站的高程差除以上地壳波速近似得到垂向的相位差。这样最终确定了各台站相对于参考台站（2号）的相位差，以此对各个台站记录的相位予以调整。

由于大地震中近场地震动峰值高且中低频成分丰富，同时含有较多的低频成分（周期1 s以上）（李小军，2001），再考虑到对于大地震的近场成像研究，采用低频成分往往能够较准确地估测震源尺度（Meng et al，2012 ），因此本文对地震动时程进行0.2—1 Hz的带通滤波，最终得到的加速度时程如图4所示。在3 s时间窗下，以1 s为滑移间隔，通过寻找每一秒的相关系数和的峰值位置来确定破裂主能量源方位角。预先定义龙门山断裂带为过震中、走向为225°的线型断层，则自贡台阵朝向该方位角的射线与这个预定义断层的交点即为破裂主能量源的位置。

图  4  经过震相对齐、滤波和相位调整后的地震动加速度时程

Figure  4.  The time histories of ground motion after seismic phase alignment，filtering and station delay processing

下载: 全尺寸图片 幻灯片

图5给出了离线条件下汶川地震的破裂长度估测结果随时间的变化情况。图5a中，颜色代表断层各位置的相关系数和，圆点的大小代表每一时刻叠加的相关系数峰值，圆点越大说明能量越集中。地震波刚刚到达时，能量的来源集中，相关系数峰值较大，而加速度时程在100 s后很少出现大的相关系数峰值，这是由于这个阶段各种传播路径的地震波的能量叠加、成分和能量来源非常复杂所造成。

图5b为参考台站（2号）的垂向加速度时程，可以看到，地震动加速度时程显示在约18 s时Pn波传播至该台阵，约45 s时Sn波到达，约78 s时地震动时程达到了其加速度峰值（其它台站记录的地震动加速度峰值均在该时刻附近），此时图5a中估测出的破裂断层从震中向东北和西南方向各传播60 km、总长约为120 km的双向破裂，得到估测结果的时长为65 s。

图  5  离线条件下汶川地震瞬时主能量源实时估测结果（a）及同时刻2号参考台站的垂向加速度时程（b）

图（a）中圆点的位置代表这一时刻的主能量源位置，圆点的大小代表这一时刻相关系数和的极值，虚线表示两阶段估测的破裂长度范围

Figure  5.  The outline result of the real-time estimation of the main energy source of Wenchuan earthquake （a） and the vertical acceleration record of the reference station No.2 at the same time （b）

In Fig.（a），the dot position represents the location of main energy source at the moment，the dot size represents the maximum size of correlation coefficient at the moment，and the dotted lines represent range of the rupture

下载: 全尺寸图片 幻灯片

从图5还可以看到，65 s时观测点的地震动远未达到该地的地震动最大值（位于70—80 s），说明在地震动加速度峰值到来之前便可获得较合适的破裂长度估测结果，这一先于地震动最大值出现之前的破裂长度估测对于地震预警系统非常重要。进一步可看到，120 km的断层破裂长度是由45—65 s时段的主能量估测结果所决定，反映了汶川大地震最主要的一次能量释放的破裂长度，该结果与张勇等（2008）的结论相似，表明本文方法能有效地给出最主要破裂过程的估测结果。根据张勇等（2008）的研究结果，汶川地震最主要的能量释放事件的发生时间介于震后14—34 s之间，该过程所释放的地震矩约占整个破裂过程地震矩的60%，相应滑移量超过3 m的断层长度约为120 km。考虑到地震波传播的时间因素，我们认为图5中45—65 s的估测结果较好地描述了发震后14—34 s这一阶段的断层破裂长度。此外，观察到该方法容易低估远离台阵方向上的破裂长度。120 s时利用本文方法得到的断层破裂长度为230 km，其中向南西延伸约170 km，向北东延伸约60 km，这个结果与其它研究结果相比有一定的差距，除了较张勇等（2008）的结果呈现较为明显的低估之外，还缺失一些震中东北方向上的重要破裂过程。从图5的结论来看，破裂过程到达42 s时，震中西南方向约150 km处出现了新的破裂，几乎在同一时间，虽然发生在北东方向上的新破裂的滑动量更大，但是这一过程未被探测出来。本文认为是震中北东方向上的破裂距离台阵较远，能量衰减较多，且该阶段中传播到达台阵的地震波成分复杂所致，所以北东方向上的破裂没有显现出来，而西南方向上的破裂距离台站较近，破裂过程被记录了下来。Meng等（2014）关于2012年El Mayor-Cucapah地震的研究中也出现了类似的现象。因此，这一方法对于汶川大地震主要破裂过程的破裂长度能够得到较好估测结果的同时，会在后期的估测过程中出现离台站较远的北东方向上破裂长度呈现明显的低估。这一现象与台阵相对于发震断层的距离较远有关。若在震中的北东、西南方向上有相应的台阵记录，将能够解决破裂长度低估的问题。这就是基于地震动观测台阵记录开展大地震断层破裂长度实时估测所面临的的限制条件。

本文所研究的断层破裂长度估测方法在地震预警原地、异地信息检测模式上有一定的局限性。原地地震预警系统是根据P波波速比更具破坏性的S波、面波波速快的物理原理，利用P波来估测未来S波、面波的影响程度。从图5中可以看到，本文的这种方法在S波到来后才能得到较为科学的估测结果，所以难以在原地地震预警系统中发挥作用。而异地地震预警系统是根据离断层较近的估测结果来为破坏性地震波还未到达的地区提供预警信息，本文方法能够基于近断层的地震动观测台阵记录计算得到断层破裂长度信息，以降低预警系统低估地震动强度和强烈地震动分布范围的风险。

如果要在S波到来前对震源长度进行有效的估测，可以从早期P波中寻找信息。如图5所示，在断层破裂长度的估测值达到120 km时，虽然地震动峰值加速度还未出现，但S波首波在约20 s前已经到达，这对于大地震来说是相当危险的，尤其是在近场，地震动峰值大，S波到达后很快会造成严重的破坏。根据早期的P波信息所得到的断层破裂长度为从震中向南西方向传播约50 km，由于这种方法对初始震中位置的估测出现了偏差，所以在破裂方向的估测上同样出现了偏差；但破裂过程得到了较为完整的记录，即从初始震中向北东方向传播，再从北东方向向南西方向传播。如果利用P波信息对汶川地震的主要破裂过程进行推断，结合初始震中的位置，可以得到汶川地震的初期破裂长度为从震中向北东方向延伸50 km，这虽然与利用S波进行估测的结果相比较小，但得到估测结果在时间上要比利用S波的早约40 s，而且这一结果也能够反映出地震破裂传播的主要方向。现阶段采用近场P波信息对震源破裂长度的估测研究还较少，这种估测方法仍需更多的大地震台阵记录进行准确性验证。

震源破裂长度的估算若要满足预警系统的需求，除了在危险来临前获取计算所需的数据，对算法本身的计算速度也有较高的要求。本文的时间窗漂移量设为1 s，那么，对于每一个时间窗的计算时间应在1 s以下，甚至要远小于1 s，以便为数据读取、通讯预留时间。

对算法本身的计算速度进行测试分析，具体的计算环境为：Intel i7-2600 k中央处理器，4核8线程，最大主频3.40 GHz；可用内存3.47 G；32位win7操作系统；Matlab2011a软件；并行算法的计算节点为8个。

算法速度的测试结果显示：数据准备过程用时4.149 7 s，利用断层实时估测的串行算法计算每1 s的断层破裂位置用时不超过2.4 s。对计算速度有显著影响的因素为：① 台阵的台站数量n_s，用于计算每两个台站间的相关系数，时间复杂度为O（n_s ²）；② 断层的分割数量n_f，时间复杂度为O（n_f）；③ 慢速矢量的间隔数量l_s，时间复杂度为O（l_s）；④ 时间窗长度n_t。总的计算量为O（n_s ²·n_f·l_s·n_t），通过对整个断层分割进行8线程的并行计算，能够有效地提高计算速度。在并行计算初期，CPU需要经过一个加速的过程，之后该算法得到每一秒估测结果的时间均能够稳定在0.8 s附近，因此该算法理论上能够满足实时断层长度估测的应用。

3.   讨论与结论

本文采用近场地震动波束形成技术，利用自贡地形地震动观测台阵在汶川地震中得到的地震动时程，通过反演破裂主能量源方位角的方法开展了汶川地震发震断层破裂长度的实时估测，主要结果如下：

1） 完善了利用近场小口径地震动观测台阵实时估测断层破裂长度的技术路线，开发了用于估测断层破裂长度的算法，提出了调整各台站相位的方法，即在每个台站的延时无法确定的情况下，仍能够近似地得到每个台站之间的相位差。

2） 在地震动峰值加速度到来前约10 s，该方法得到的断层破裂长度为120 km，从震中向北东、南西各传播60 km，这一结论较好地反映了汶川大地震最主要的能量释放过程的断层破裂长度。

3） 受能量衰减因素的影响，该方法容易低估远离台阵方向上的破裂长度。最终的破裂长度估测结果约为230 km，主要从震中向南西方向延伸，这一结果低估了北东方向上的破裂，本文认为可通过在震中的北东、南西方向上各放置一个台阵的办法来降低破裂长度低估的可能性。

4） 讨论了使用近场P波估测破裂过程的可行性，通过该方法得到汶川地震的初期破裂为从震中向北东方向延伸50 km，且得到估测结果的时间比利用S波的早约40 s。基于P波信息的方法能够较准确地描述断层的初期破裂传播过程，但是对震中的估测存在一定的偏差。

5） 为提高实时估测震源破裂过程方法的实用性，分析了算法的计算速度。采用8线程并行计算算法，该算法下计算每秒的破裂位置用时约为0.8 s，能够为数据传输预留时间并满足实时破裂长度估测的应用。

本文的研究进一步表明，破裂长度估测方法仍然有其局限性，因此存在改进之处：第一，考虑到本文所利用的地震动观测台阵获得的强地震动记录数据较少，本文在处理台站之间高程带来的相位差时采取了近似的方法而不是利用大量数据进行拟合分析，这会给加速度时程的相位确定带来偏差；第二，相比美国的UPSAR台阵，自贡台阵的观测点数量少，台阵观测点口径也较小，且自贡台阵与龙门山断裂带相距较远，因此在计算主能量源位置时，空间不确定性较大，容易引起对破裂位置估测的偏移；第三，该方法虽然能够将主要的破裂过程估测出来，但是当一些局部破裂发生在远离台阵的断层走向上时，不容易被探测到，这会导致断层长度的低估。尽管如此，该方法仍然能够在大地震发生后第一时间内为其破裂长度的估测提供重要信息。