2013, Vol. 56
2. 中国地震局地震预测研究所, 北京 100036;
3. 天津市地震局, 天津 300201;
4. 中国科学院地球深部研究重点实验室, 中国科学院地质与地球物理研究所, 北京 100029
2. Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China;
3. Earthquake Administration of Tianjin Municipality, Tianjin 300201, China;
4. Key Laboratory of the Earth's Deep Interior, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
我国历史资料中有着大量历史地震破坏程度的记载.与7级以上的大震或巨震相比，对6级左右中等强度历史地震的震源参数研究相对并不充分.利用历史资料定量估计得到较可靠的6级左右历史地震震源参数，可有效弥补现代仪器记录的地震资料相对短暂之不足，为强震危险性分析提供重要的基础资料.
由历史资料记载的地震宏观破坏记录确定相应烈度分布，并依据高烈度区的几何中心、极震区和有感区面积及烈度分布形状来估计历史地震的震源参数（多为震中经纬度），是普遍使用的方法[1-2].但因当时人文、经济和地理环境差异导致史料记载的震害详简程度不一，由残缺不全的历史记录推断的6级左右地震烈度数据点相对稀少，且空间分布图像不大完整，由此勾画的烈度等震线分布不免带入研究者的主观假设，所给出的震中参数存在相当大的主观性和不确定性.
为克服利用烈度数据估计震中位置和震级主观性强以及不确定性无法定量估计的困难，Bakun和Wentworth[3]以美国加州地区为例，用已知震中位置和震级的现代仪器记录地震作为训练集，确定震级-震中距-烈度值衰减关系式，使用格点搜索试算法估计震中和震级，并通过对训练集数据反复抽样给出震中和震级估计的置信区间.Bakun等[4]近来又对方法做了改进，不通过训练集而直接对地震烈度数据重采样，使用Bootstrap方法估计震中和震级的置信区间.Bakun-Wentworth方法直接对烈度数据点拟合而放弃烈度分布图，具有一定的可重复性和客观性，因而得到了广泛的应用，但其局限性在于必须依靠一定数量的现代中强地震来确定衰减关系，对如我国华北那样历史地震资料丰富，而有现代仪器记录之破坏性地震相对缺乏的地区，Bakun-Wentworth方法应用较为困难.
通过理论计算不同震源参数在地表产生的烈度分布，并与烈度等震线图相比较，也是估计历史地震震源参数的一个途径.陈培善[5]早在20世纪70年代就对利用位移表示定理的中远场近似计算近场值进行了探索.Sirovich等[6]通过枚举可能取值空间的设定震源参数，利用S波远场辐射花样计算各种设定震源产生的烈度分布图，与历史资料比较确定震源参数的方法（称为KF方法），具有可估计多个震源参数，且不依赖现代中强地震记录的优势.但KF方法仅使用S波远场辐射花样，没有考虑地震波中体波低频部分和面波的影响，计算出的双峰型烈度等震线与一般高烈度区近椭圆烈度分布形态符合程度不佳，并且在拟合Ⅵ度以下的烈度等震线时情况较差[6].
随着计算机技术的不断发展，利用地震波场模拟计算估计历史地震震源参数的应用研究越来越多[7-10].本研究将Bakun-Wentworth方法直接拟合烈度数据点和KF方法参数枚举的优势结合起来，引入近年来在震源破裂过程反演[11]和地震危险性分析[12]得到应用的模型选择（model selection）和多模型统计推断（multimodel inference）方法[13]，发展一种利用地震波场数值模拟和烈度数据估计多个震源参数及其不确定性的方法，克服KF方法仅使用S波高频近似的问题和Bakun-Wentworth方法难以在现代中强地震活动记录较少地区的使用问题.
2 模型选择方法和算法流程本方法将所有可能的震源参数组合进行枚举，每个组合视作一个震源模型，通过模型选择方法对所有可能的震源模型做出评估，给出各震源参数的权重值分布.
2.1 模型选择的数学基础模型选择方法通过计算各个模型与观测数据之间的相对熵（relative entropy），对模型拟合数据的程度给出定量化估计.这里简要介绍本研究所需的统计量及其表达式，对模型选择理论详细系统的论述见文献[13].
首先计算出观测或推断的数据点烈度值与理论模拟计算得到的烈度值间的方差
|
(1) |
其中ISi和IBi分别为第i个点的模拟计算烈度值和烈度记录值，N为烈度数据点总数.接着采用适用于小样本情况的赤池信息准则[13-14]来计算每个设定震源模型的AICc值，即各震源模型与烈度记录间的相对熵
|
(2) |
其中K为模型参数个数，例如本研究考虑震中经度、纬度、震源深度和滑动角的情况下，K=4.在所有设定震源模型的AICc值中找出最小的AICcmin，进而计算在给定某些参数下模型子空间中模型的AICc差值Δi
|
(3) |
根据AICc差值可得到相应模型权重wi
|
(4) |
其中构造底数La=1.5.如此得到的wi在相应震源模型参数空间的和为1，且每一个wi都不小于零.因而wi满足完备性、非负性和可列可加性，由概率分布的公理化定义可将wi理解为一组概率分布.由此可对具体参数进行多模型统计推断，即依据模型权重计算参数的权重值分布，估计相应的震源参数.其计算方法类似于对多维离散随机变量的全概率做边缘化（marginalization），求某一参数维度的边缘概率分布.
2.2 方法计算流程本文考虑到发震断层走向和倾角可充分借鉴地质构造调查、人工地震探测和烈度几何分布等先验估计，故仅关注震中位置、震源深度和滑动角构成的四维震源参数模型空间.所采用方法首先计算震中位置搜索格点AICc值分布，并据此选择可能的震中位置，再根据每个已选震中位置估计震源深度和滑动角.方法计算流程如下：
①震源深度、滑动角和发震断层破裂方式在可能取值范围以一定间隔进行枚举，组合成为备选震源模型.
②计算每个备选震源模型在地表的理论地震波场，依据水平向峰值加速度（PGA）或峰值速度（PGV）与烈度的换算关系（如表 1所列的中国地震烈度表GB/T 17742-2008或相应研究区的关系式）得到烈度分布.考虑到本研究使用的时域伪谱法（PSTD）[15]理论地震图计算直接给出的是速度记录及速度与地震波能量的相关性[7]，故本文采用将PGV转化为烈度值的方式，而不采用将理论计算的速度记录转化为加速度再利用经验关系转为烈度值的方式.
|
表 1 中国地震烈度表（GB/T 17742-2008） Table 1 China earthquake intensity scale given by GB/T 17742-2008 |
③首先利用式（1）和（2）计算每个震源模型的烈度方差及其AICc值，然后取每个震中位置相应模型的最小AICc值作为此震中格点的AICc值，得到震中搜索范围内的AICc值分布.
④根据震中AICc值分布，选择AICc值较小的格点作为震中位置模型选择初步结果，并采用式（3）和（4）计算震中格点的权重值.
⑤在步骤③得到的震源模型AICc值基础上，对初步选择的每个震中位置结果，依据相应的AICc值利用式（4）分别求得震源深度和滑动角的权重值分布，通过多模型统计推断震源模型的震源深度和滑动角估计结果.
⑥结合地质构造、地震学和大地测量等研究成果，对以上震源参数模型选择结果的合理性进行分析取舍，最后得到历史地震震源参数估计的可能结果.
3 方法数值检测为检测方法复原已知震源参数的能力，本研究基于表 2所示的一维速度模型，将已知的震源参数作为双力偶点源，采用时域伪谱法（PSTD）[15]来计算理论地震图，并将理论波形的水平峰值速度（PGV）转化成烈度值.所使用的转化关系参考中国地震烈度表（表 1），利用中值回归得到烈度值与PGV的关系式为I=3.26logPGV+3.45，其中PGV单位是cm/s.对转换所得的烈度值加入一定标准差σ生成的随机噪声扰动值后（取四舍五入后的整数），作为烈度数据点来估计震源参数.
|
|
表 2 数值检测使用的一维速度结构模型 Table 2 1-D structure used in this numerical test |
鉴于我国6级左右历史地震的烈度数据点记录多为几十个，故不失一般性以生成随机分布的30个烈度数据点来进行数值检测.所检测震源参数的震中位置搜索格点间隔为2 km，震源深度在3~20 km范围以1 km为间隔取值，滑动角在0°~350°范围以10°为间隔取值（滑动角取值规则参见文献[16]）.
对简单的一维速度结构，所检测震源参数模板的理论震中可取在正中间的位置，而理论震源深度和滑动角则分别参照1998年张北6.2级地震和2006年文安5.1级地震震源参数来设定，即分别取深度10 km和滑动角130°[17]与深度9 km和滑动角210°[18].将参照的文安地震震源深度10 km改为深度9 km，是为了检测不同震源深度模板的情况.这两个设定的震源参数模板，一个是以逆冲为主兼有右旋走滑，一个以右旋走滑为主兼有正断分量.检测中在理论计算转换的烈度值中加入均值为0、标准差σ分别为0.5和1的高斯分布随机数作为观测值，来模拟因地下三维速度结构和场地效应及人为因素等导致的烈度数据的不确定性.
Bakun等[4]利用Bootstrap重采样1000次计算分析历史地震震源参数估计方法的稳定性和误差.本研究同样对两组震源参数模板分别进行1000次参数估算，图 1和图 2分别是两组模板参数测试的统计分析结果.其中图 1a和图 2a以直方图形式给出了理论转换烈度加入不同标准差的随机干扰后的烈度值与理论值的偏离情况，可见考虑了烈度偏差达到1度甚至2~3度的情形.图 1b和图 2b给出了以最大权重作为震中估计结果的频次叠加图，图 1c和图 2c则给出了以最大权重估计的震中格点与理论震中的偏离统计（统计间隔为震中搜索间隔的2 km）；由此二图可见重复多次以最大权重估计的震中结果都聚集在理论中心点周围，与理论中心点的距离越大（大于4 km）而被作为估值的频率越低.考虑到历史地震资料和震中估算精度的不确定性，图 1d和图 2d统计给出了以权重值超过30%最大权重的格点作为震中估算结果的频次分布，即表示格点在1000次模型选择测试中被选中的次数；由图可见，即使在扰动标准差σ=1时，理论震中位置格点被选中的频次分别达912次和920次，虽然可选的震中格点增多且聚集范围有所扩大.1000次随机生成烈度点并添加一定干扰值的震中测试，模型选择的结果大量聚集在理论震中位置及其周围，说明本文所采用的分析方法对震中位置估计在统计上具有一致性和统计无偏性，以5 km估算精度统计的震中选择结果的可信度超过了90%.
|
图 1 震源深度10 km和滑动角130°的震源模板1000次测试的统计分析结果 （a）加入不同标准差的随机扰动后的烈度偏差统计直方图；（b）只选取最大权重值推断震中位置格点的频次叠加图；（c）震中估计结果中最大权重格点与理论震中位置的距离统计；（d）取权重值大于最大权重的30%的格点作为震中估计结果的频次分布图；（e）震源深度的最大权重估计值与理论值的偏差统计柱状图及其权重值分布；（f）滑动角最大权重估计值与理论值的偏差统计柱状图及其权重值分布. Fig. 1 Statistical results of the numerical test with input parameters of focal depth 10 km and rake 130° (a) Statistical distribution of intensity shift; (b) Statistical results of epicenter dislocation distance between maximum weight grid and the input location; (c) and (d) are frequency distributions of epicenter grids from model selection by the maximum weight and large than 30% of maximum weigh trespectively; (e) and (f) are weight distributions and misfit histograms between maximum weight solution and input value of focal depth and input value of rake angle respectively. |
|
图 2 参数为震源深度9 km和滑动角210°的模板1000次测试统计分析结果 （a）-（f）同图 1. Fig. 2 Statistical results of the numerical test with input parameters of focal depth 9 km and rake 210° Others are sameas Fig. 1. |
图 1e和图 1f与图 2e和图 2f分别是震源深度和滑动角的相应统计结果，中间的累积柱状图分别表示以最大权重估计的震源深度和滑动角与理论值的偏差数，不同颜色表示不同的偏差值；左右二图展示的震源深度和滑动角权重值分布显示：理论震源深度和理论滑动角的估计值权重最大，说明本文所采用的分析方法对震源深度和滑动角的估计在统计上也具有一致性和一定的无偏性.
4 方法检测实例：2004年美国Parkfield6.0级地震在数值检测结果的基础上，我们进一步以获得较多强震记录的2004年美国加州Parkfield6.0级地震为例来检测该分析方法的可行性.2004年Parkfield6.0级地震发生在右旋走滑的San Andreas断层Parkfield段，地震破裂长度约为40 km，由起始破裂点向北西和南东方向不对称双向破裂，向北西破裂约30 km，向南东破裂约10 km，破裂速度约3 km/s[19-20].Berkeley大学的快速矩心矩张量解（http：//quake.geo.berkeley.edu/mt/nc51147892.mt）给出的该地震震源深度为8.0 km和滑动角为185°（-175°）.全球GCMT解（http：// www.globalcmt.org）给出震源深度12 km和滑动角188°（-178°）.利用强震波形和GPS数据联合反演震源破裂过程[19, 21]得出的震源深度8.1 km.
我们从COSMOS强震数据中心（http：//db. cosmos-eq.org）收集到如图 3a和表 3所示的Parkfield地震的39个强震台站记录波形，并将强震波形记录的水平向峰值速度采用表 4的对应关系转化为烈度值，检测本文方法的参数估计结果与利用波形反演得到的震源参数是否一致.
|
图 3 Parkfield地震强震观测点与震源搜索模型示意图 （a）黑色三角为强震台站位置，也是本文使用的烈度数据点位置，虚线网格为震中搜索格点.（b）三种破裂方式示意图，分别是对称双向破裂和向两侧的不对称双向破裂. Fig. 3 Intensity data distribution and epicenter search grid of the Parkfield earthquake (a) Black triangles are strong motion stations and also intensity data points locations used in the test, epicenter search grids are given by dash lines.(b) Three kinds of rupture styles with symmetric and asymmetric bilateral rupture in two directions. |
|
|
表 3 本文使用的强震台站名称及经纬度 Table 3 Strong motion stations used in the Parkfield earthquake test |
|
|
表 4 PGV转化为烈度值对应表1） Table 4 Relationship between instrumental intensity and PGV from USGS1) |
本研究的理论地震图计算以San Andreas断层为界采用不同的一维地下速度结构模型（表 5），即已用于Parkfield地震震源破裂过程反演[19, 21]的速度模型.
|
|
表 5 Parkfield地震地震动模拟计算所用的一维速度结构模型 Table 5 Crustal velocity structure used in the Parkfield earthquake test |
我们同样采用Ma等[19]的发震断层垂直于地表和断层走向140°的简化假设，设断层破裂长度为40 km，并假设在整个断层面上地震滑动角一致，断层面上应力降分布为从起始破裂点向周围衰减.备选震源模型构建考虑三个因素：震源深度、滑动角和断层破裂方式.震源深度取值范围3 km到15 km，以1 km为间隔枚举，滑动角取值范围0°到350°，以10°为间隔枚举.如图 3b断层破裂方式分为对称双向破裂和分别向两侧传播的不对称双向破裂，起始破裂点分别取在断层北西端的1/4、1/2和3/4处，称为L方式、M方式和R方式.图 3a中网格为震中估计格点搜索范围示意图，搜索格点间隔为5 km，震源球所指位置是GCMT解给出的震中位置.
图 4给出了Parkfield地震震源参数估计过程的主要计算结果，其中图 4a为震中搜索格点的AICc值分布；图 4b为震中模型权重值分布，权重最大的震中格点为图中模型1所示的震中位置，其权重值为0.367；以不低于最大权重值的30%来选择，则有如图所示的4个震中位置估计模型，即模型1至模型4（权重值相应为0.367、0.16、0.108和0.108）.其中模型1权重值明显大于其它三个，且与图 3a中GCMT解给出的震中位置一致.
|
图 4 Parkfield地震震源参数估计结果 （a）震中格点搜索AICc值分布图；（b）震中模型选择结果，圆圈表示选中的4个模型震中位置；（c）以模型1为震中估计震源深度和滑动角的AICc值分布图；（d）震源深度的权重值分布；（e）滑动角的权重值分布. Fig. 4 Source parameters estimation result of the Parkfield earthquake (a) AICc distribution of epicenter search grids; (b) Epicenter model selection results with 4 circles; (c) AICcdistribution with different focal depths and rake angles for the epicenter Model 1;(d) and (e) are weight distributions with focal depth and rake angle, respectively. |
以模型1的震中位置为依据，进一步估计Parkfield地震震源深度和滑动角.图 4c是震中位置取模型1位置时，不同震源深度和滑动角模型的AICc值分布，图 4d和图 4e分别是震源深度和滑动角的权重值分布.结果显示震源深度在7 km处取得最大权重值0.29，8 km深度权重值略小为0.27，这一结果与Berkeley快速矩心矩张量解震源深度8.0 km和震源破裂过程反演[19, 21]得出的8.1 km震源深度结果基本一致，但与GCMT解给出的12 km深度结果有差距.滑动角权重值分布曲线有两个峰值区，一个代表的震源机制是左旋走滑兼有逆冲分量，另一个的震源机制则是右旋走滑兼有正断分量，需要根据地质构造等其它信息选择其中一个.考虑到San Andreas的右旋走滑断层属性，及1966年同一断层段发生的6级地震的震源机制为右旋走滑并带有较小的倾滑分量[22-23]，因此可以判定2004年Parkfield地震震源机制应取右旋走滑兼有正断分量的峰值，含左旋走滑分量的峰值被舍去.滑动角在200°取得最大权重值0.23（舍去一个峰值则权重值比图 4e上的读数增加一倍），在170°、190°和210°权重值相近约为0.11，220°和230°权重值也超过了0.05，因而滑动角估计结果为[170°，230°]，震源机制以右旋走滑为主，兼有少量正断分量.此结果与Berkeley快速矩心矩张量解滑动角185°（-175°）和GCMT解滑动角188°（-178°）的结果基本一致.
取模型选择结果中权重值最大的震中位置、震源深度8 km和滑动角200°，计算了39个台站的理论地震图（图 5）与实际记录波形进行对比，两个水平方向分别为平行断层和正交断层方向.波形长度取15秒，经过0.2~1 Hz的4阶零相位Butterworth滤波器带通滤波.模拟波形与实际记录波形的符合程度是可以接受的，说明震源模型选择结果与实际发震断层性质相近.
|
图 5 Parkfield地震水平向模拟波形和记录波形对比图 红色为模拟波形，黑色为记录波形.两个方向分别为平行断层和正交断层，波形经过0.2~1 Hz的4阶零相位Butterworth滤波器带通滤波. Fig. 5 Waveforms of horizontal components from observed ground motions (black) and synthetic calculations (red) with the best selected source model Left and right seismograms of each station are components parallel to and perpendicular to the fault.All the waveforms are band-pass filtered between 0.2 to 1 Hz with a zero-phase forth order Butterworth filter. |
需要说明的是，由于本研究是以Parkfield6.0级地震进行应用检测，震中位置已有不少研究给出较确切的结果，故仅以模型1的断层位置进一步估计震源深度和滑动角.而在应用于具体的6级左右历史地震震源参数估计时，则需对出现的多个震中位置估计结果，分别推断震源深度和滑动角，得到多组可能的估值.最后根据区域地质构造和大地测量等研究结果，综合判断这几组估值的合理性.
5 方法应用实例：1882年河北深县6级地震本研究进一步以有着较好历史地震破坏记录的1882年河北深县6级地震为例，采用同样的途径对深县地震的震源参数进行估计.《中国历史强震目录》[1]中给出的深县地震（图 6a）震中位于38.1°N、115.5°E（震中精度属于2类，误差≤25 km），极震区在深县北部，较高烈度的Ⅷ度区和Ⅶ度区位于深县凹陷区内.梁富康等[24]据反射勘探剖面编制的深县凹陷区断裂体系分布图（图 6b）为我们的应用研究提供了构造参考依据.梁富康等[24]将深县凹陷区的断层规模分成三级，即2条一级断层（衡水断层和虎北断层），4条二级断层，从北到南依次为旧城北断层、何庄断层、深西断层和深南断层，及广泛发育的三级或级别更小的断层.梁富康等[24]提供的通过深县凹陷主体的构造剖面显示，深县凹陷区主要断层为正断性质，区内沉积层厚度在5 km左右.衡水断层为走向北西、倾向北东、具有左旋走滑分量的正断层，断层倾角较陡，一般55°左右，有些地段断面陡直[24-25].旧城北断裂走向近东西，倾向南，何庄断层也为东西走向，向南倾斜[24].虎北断层为北东东走向、北北西倾向，深西断层为北东东走向、南南东倾向，深南断层也是北东东走向，而其倾向为北北西[24].
|
图 6 《中国历史强震目录》给出的1882年深县6级地震烈度等震线图（a）和深县凹陷区断裂体系分布图（b，改自梁富康等[24]） Fig. 6 Isoseismal map of the Shenxian earthquake given by China historical earthquake catalog (a) and regional faults systems in the Shenxian depression (b, modified from Liang et al.[24]). |
龙峰等[26]以系统收集整理获得的1965年以来华北地区发生的34次中强地震的可靠破裂尺度参数，回归得到华北地区6级左右地震断层破裂长度在15 km左右.据此判断，图 6中的三级以下断层规模无法成为深县6级地震的发震断层，4条二级断层和2条一级断层有可能是发生深县6级地震的断裂.据深县地震的烈度等震线形态判断，发震构造似乎应为北西-南东走向的断层，但北西向的衡水断层与极震区距离在20 km以上，且处于Ⅶ度区以外，似乎不具备成为深县地震发震断层的条件.因而1882年深县地震的发震构造是一个值得研究的问题.
5.1 烈度资料依据前人整理的《河北省地震资料汇编》[27]和《中国历史强震目录》[1]，深县地震给出有破坏程度记载的地点名称42个，深县、深泽等较大城镇在烈度图（图 6a）上有明确标识，可以以此为依据确定烈度值；高烈度区及周边的村庄则根据周围城镇的烈度值和文献记载中关于“轻”和“重”的描述来判定烈度值.深县处于Ⅷ度和Ⅶ度区之间，则其中“最重”和“较重”的村判为Ⅷ度和Ⅷ度弱（郭家庄、北安、深州城关四厢、西八弓），“较轻”和“又轻”的村判为Ⅶ度（西马庄、南官村、北官村）.束鹿处于Ⅶ度和Ⅵ度区之间，其中“最重”和“次重”的村判为Ⅶ度（双井、北庞营、南庞营、北吕彩），“较轻”的判为Ⅵ度（沙河）.安平的烈度为Ⅵ+，郭家店（郭屯）是周围村庄中较重的，因此判为Ⅶ度.对于文字资料没有单独明确描述的，参考烈度等震线图判断烈度.“城北之郭家庄等村，城西之杜家庄等村，西北与安平、束鹿毗连之张村等村，同时震塌房屋，间有伤毙人口等情”.根据烈度等震线图，张村在深县西北，与图中给出的Ⅷ度区很接近，与郭家庄一致判为Ⅷ度，杜家庄在深县正西，与Ⅷ度区相距较远，判为Ⅶ度.破坏程度记录不明确和无法查明确切经纬度的烈度点舍去.最终使用的深县地震39个烈度数据点位置及烈度值见图 7和表 6.
|
图 7 深县地震烈度数据点分布图 Fig. 7 Intensity data points distribution of the Shenxian earthquake |
|
|
表 6 1882年深县地震烈度数据点位置分布 Table 6 Intensity data points of the 1882 Shenxian earthquake |
根据深县地震的烈度数据点及其等震线和区域构造分布（图 7和图 6），震源参数估计以三个可能的断层走向来分析，分别为与衡水断层相近的北西走向（走向设为135°）、与旧城北断层和何庄断层相近的近东西走向（设为90°）、及与深南断层、深西断层和虎北断层相近的北东东走向（设为60°）.理论地震图计算采用表 2所示的速度模型，参考中国地震烈度表（表 1）将理论波形水平峰值速度（PGV）转化成烈度值，作为烈度数据点来估计震源参数.
震源参数模型的地震震中搜索围绕39个烈度数据点中心位置（115.677°E、38.179°N）的100 km×100 km范围进行，搜索间隔5 km，格点的排列方向与可能的发震断层走向一致；震源深度则在3~20 km之间以1 km为间隔枚举，断层倾角在30°~90°之间以15°为间隔枚举，滑动角在0°~360°间以10°间隔枚举；断层破裂长度设为20 km，破裂方式分为从断层两侧开始的不对称双向破裂和从断层中间开始的对称双向破裂三种破裂模式；考虑到历史地震震级估计的不确定性，以及浅部沉积层对地震波的放大作用，震级在5.5~6.5级间以0.05级为间隔枚举.
图 8-图 10分别给出了以北西走向、东西走向和北东东走向断层为可能的发震断层的震源参数估计结果.发震断层为北西走向的震中格点搜索结果（图 8）中，有12个格点的AICc为最小值-77.0，另有7个格点AICc值为次小的-68.3，这19个格点作为震中位置模型选择的初步结果，如图 8b中圆圈位置.考虑到深县凹陷区的沉积层厚度及地震较少发生在沉积层内，因此舍去对应震源深度最大权重值小于5 km的模型.若不同的震中格点对应模型为相同倾角和破裂模式，则只选择其中一个AICc值最小的格点为此倾角和破裂模式对应的震中位置.则有4个格点成为北西走向发震断层震中模型选择的最终结果，称作模型1-模型4，模型2和模型3震源机制为正断为主兼有走滑分量，模型1和模型4以左旋走滑为主.
|
图 8 北西走向断层震源参数估计结果 （a）震中格点搜索AICc值分布；（b）震中模型权重值分布和震中模型选择初步结果，图中圆圈是初步选择的震中位置；（c）震中模型选择的可能结果，不同颜色标出了可能选择模型的震中位置与发震断层位置和破裂情况，对不同模型有所重叠的断层稍作偏移以清晰展示结果；（d）-（g）相应震中选择模型的震源深度和滑动角估计结果. Fig. 8 Source parameters estimation results with NW direction strike focal structure (a) AICc distribution of epicenter searching grid; (b) Model selection result of preliminary epicenters, the circles demonstrate preliminary model selection epicenter locations; (c) Model selection result of preferred epicenters named Model 1-4 with relative fault and its rupture style; (d)-(g) Multimodel inference results of focal depth and rake angle for four models. |
|
图 9 东西走向断层震源参数估计结果 （a）-（g）意义同图 8 Fig. 9 Source parameters estimation results with EW direction strike focal structure Others are same as Fig. 8. |
|
图 10 北东东走向断层震源参数估计结果 （a）-（k）意义同图 8. Fig. 10 Source parameters estimation results with NEE direction strike focal structure Others are same as Fig. 8. |
对东西走向发震断层设定的震中格点搜索AICc值分布（图 9a），共有6个格点AICc取得最小值-68.3，13个格点取得AICc次小值-61.2，这19个格点作为震中位置模型选择的初步结果，如格点权重值分布图 9b中所示圆圈位置.对初步选择的每个震中位置，使用多模型统计推断估计震源深度和滑动角（图 9d-9g）.同样对震中初步结果做多模型统计推断，得到4个震中位置作为东西走向发震断层震中模型选择的最终结果（图 9c），称作模型5-模型8，震中位置比较集中于深县以北15 km范围内.考虑到研究区内主要断层的正断属性，因而在滑动角两个峰值中取含有正断分量的部分（180°~360°），4个模型都以正断为主兼有少量走滑分量.
发震断层为北东东走向震中位置估计初步结果如图 10a，有2个格点的AICc为最小值-77.0，另有9个格点AICc值为次小的-68.3，这11个格点成为震中位置模型选择的初步结果，5个模型成为最终的模型选择结果，称作模型9-模型13.震中位置分布比较集中，其中模型10、模型11、模型12和模型13震源机制以正断为主，模型9兼有正断和左旋走滑分量.
5.3 模型选择结果讨论
表 7所列的模型选择得到的13个结果中，烈度误差和表达式为
|
|
表 7 深县地震模型选择结果汇总 Table 7 Model selection results of the Shenxian earthquake |
图 11a为断层走向北西的模型震中和断层位置，模型1和模型4震源机制估计结果为走滑兼有少量正断分量，与区域地质构造[24-25]北西走向衡水断层正断层属性不符.模型2和模型3与衡水断裂地表出露位置相距较远，且其震源深度估计值不大于10 km，因而很难将模型与衡水断层联系起来.
|
图 11 深县地震模型选择结果综合分析图 虚线圆圈表示震中模型，内部数字表示模型号码.深县凹陷区一级二级断层位置、烈度数据点位置和三个走向发震断层震中模型选择结果等表示同图 8至图 10. Fig. 11 Source parameters model selection results of the Shenxian earthquake Dash circle demonstrate epicenter location with model ID. Others are same as Figs. 8 to 10. |
图 11b表示断层走向为东西向的四个模型震中与断层和烈度数据点的相对位置，四个模型震源机制估计结果都为正断型，与深县凹陷区内区域地质构造环境相符.其中模型7和模型8震中位置与旧城北断裂和何庄断裂位置接近，模型5和模型6震中位置在旧城北断层以北，而两条断层倾向为南南东，因而模型5和模型6震中位置似乎与东西向断层不符.
图 11c中发震断层为北东东走向的5个可能模型，震中位置都在三条断层以北.其中模型11和模型12距离深西断层较近，震源机制估计结果为正断型，也符合区域地质构造环境.
综合模型选择结果和地质构造资料分析认为：深县地震最可能的震源参数估计有两组：一是围绕模型3和模型4震中位置（震中格点中心分别为115.55°E、38.02°N，115.61°E、38.02°N），走向东西向，震源深度较深（最大权重值深度为15 km和16 km），震源机制为以正断为主；另一组是模型11和模型12震中位置（震中格点中心分别为115.51°E、38.01°N，115.57°E、38.04°N），走向北东东，震源深度较深（最大权重值深度为14 km和17 km），震源机制为以正断为主.东西向的旧城北断层和何庄断层，以及北东东走向的深西断层为1882年深县地震发震构造的可能性较大.这一结果不同于图 6a所示的主要以Ⅵ度烈度等震线形态估计的北西走向发震构造可能，其原因很可能在于对地处石家庄-保定间的定县这一烈度点值的判定和Ⅵ度等震线的勾画（图 7）.对因局部构造或场地因素或人为因素出现的单一烈度异常点，若将其视为强震动的平均衰减结果则可能导致烈度等震线图勾画的偏差.
6 结论和讨论本文借鉴Bakun-Wentworth方法[3]直接拟合烈度数据点和KF方法[6]利用参数枚举设定震源的做法，提出了一种基于烈度资料和地震动波场模拟估计历史地震震源参数的模型选择方法.数值恢复测试与Parkfield地震实例测试表明，本文提出的分析方法所估计的震源参数具有统计一致性和一定的无偏性.该方法可通过计算AICc值和模型权重，根据历史地震烈度数据估计震中位置、震源深度和滑动角等震源参数及其不确定性，克服了Bakun-Wentworth方法需要利用研究区域的现代地震资料确定衰减关系式的问题，及KF方法在拟合Ⅵ度以下烈度等震线情况较差的问题，使方法具有一定客观性和广泛的适用性.
所提出方法一个显著特点在于将利用历史地震烈度资料估计震源参数的不确定性表示为多种可能结果的形式.与现代地震的波形记录和烈度调查资料相比，由6级左右历史地震破坏资料推断得到的烈度数据不确定性较强，所含信息量较少，因而利用历史地震烈度资料估计震源参数的结果也具有较大的误差.KF方法[6]计算设定震源模拟烈度场与烈度分布图的方差，仅将取得最小方差的震源参数范围作为估计值，忽略比最小方差稍大的震源参数估值则可能抛弃潜在的结果.本文方法在分析震中格点搜索结果时，不仅选择AICc值最小的震中格点，同时也将AICc值比最小值稍大的格点也纳入震中模型选择初步结果，同时通过计算每个震中格点的权重值，尝试对震中模型的可能性做出定量化估计.将震源参数估计结果以多个可能的方式给出，充分考虑到烈度数据的不确定性，一定程度上增加了结果的可信度.
将本文提出方法应用于1882年河北深县地震的震源参数估计，基于历史资料推断的地震烈度数据点，分析认为深县凹陷区的北西向、东西向和北东东走向断层对烈度数据点拟合程度接近，三个走向的断层都有可能成为发震断层.根据模型选择结果和地质资料的综合分析，东西向的旧城北断层或何庄断层，及北东东走向的深西断层成为发震构造的可能性较大.
不可否认，受计算条件和研究程度的限制，本研究中将震源模型的发震断层设定在矩形平面，与实际地震的发震断层形态和可能较为复杂的断层破裂方式有所差距；地震动模拟计算高频地震波辐射还存在一定困难，可能会使由理论地震波场得到的烈度分布，与近断层的烈度资料相差较大；以及由历史记载的地震破坏程度来估计烈度和PGV或PGA与烈度转换的不确定性，对本文提出的方法会产生一定的影响.地下结构的横向不均匀性，尤其是浅部沉积构造导致的地震动效应，亦会影响模型选择的结果.随着地下速度结构探测研究的深入，使用高精度的三维结构模型进行波场模拟，或采用场地修正方法（sitecorrection）对烈度数据点进行合理的修正，应可提高本文方法的应用效果.
致谢感谢刘澜波教授提供PSTD理论地震图计算程序，感谢与王伟君副研究员和李乐博士等有益的讨论，感谢俞言祥研究员以及两位匿名审稿专家对本文修改完善的建议.感谢中国地震局地震预测研究所为本研究工作提供了高性能集群计算系统的支持.
| [1] | 国家地震局震害防御司. 中国历史强震目录(公元前23世纪-公元1911年). 北京: 地震出版社, 1911 . Disaster Prevention Department of China Earthquake Administration. China Historical Earthquake Catalog (in Chinese). Beijing: Seismological Press, 1911 . |
| [2] | Musson R M W, Jeménes M J. Macroseismic estimation of earthquake parameters. Sixth Framework Programme Report, 2008. |
| [3] | Bakun W H, Wentworth C M. Estimating earthquake location and magnitude from seismic intensity data. Bull. Seismol. Soc. Am. , 1997, 89(2): 557. |
| [4] | Bakun W H, Cópera A G, Stucchi M. Epistemic uncertainty in the location and magnitude of earthquakes in Italy from macroseismic data. Bull. Seismol. Soc. Am. , 2011, 101(6): 2712-2725. DOI:10.1785/0120110118 |
| [5] | 陈培善. 震源机制与烈度分布的关系(二). 地球物理学报 , 1977, 20(1): 9–19. Chen P S. Relation between the seismic source mechanism and the intensity distribution (continued). Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese) , 1977, 20(1): 9-19. |
| [6] | Sirovich L. A simple algorithm for tracing synthetic isoseismals. Bull. Seismol. Soc. Am. , 1996, 86(4): 1019-1027. |
| [7] | Fitzko F, Suhadolc P, Audia A, et al. Constraints on the location and mechanism of the 1511 Western-Slovenia earthquake from active tectonics and modeling of macroseismic data. Tectonophysics , 2005, 404(1-2): 77-90. DOI:10.1016/j.tecto.2005.05.003 |
| [8] | Nunziata C, Sacco C, Panza G F. Modeling of ground motion at Napoli for the 1688 scenario earthquake. Pure. Appl. Geophys. , 2010. DOI:10.1007/s00024-010-0113-1 |
| [9] | 周青云, 周仕勇, 陈棋福, 等. 正演推算1730年北京西郊地震的发震断层和滑动角. 地震研究 , 2008, 31(4): 369–376. Zhou Q Y, Zhou S Y, Chen Q F, et al. Forward inference of the causative fault and the slip angles of the faults of the 1730 western Beijing earthquake. J. Seism. Res. (in Chinese) , 2008, 31(4): 369-376. |
| [10] | Ding Z, Romanelli F, Chen Y T, et al. Realistic modeling of seismic wave ground motion in Beijing city. Pure. Appl. Geophys. , 2004, 161(5-6): 1093-1106. DOI:10.1007/s00024-003-2498-6 |
| [11] | Emolo A, Zollo A. Kinematic source parameters for the 1989 Loma Prieta earthquake from the nonlinear inversion of accelerograms. Bull. Seismol. Soc. Am. , 2005, 95(3): 981-994. DOI:10.1785/0120030193 |
| [12] | Scherbaum F, Delavaud E, Riggelsen C. Model selection in seismic hazard analysis: An information-theoretic perspective. Bull. Seismol. Soc. Am. , 2009, 99(6): 3234-3247. DOI:10.1785/0120080347 |
| [13] | Burnham K P, Anderson D R. Model Selection and Multi-model Inference, A Practical Information-Theoretic Approach. (2nd ed). New York: Springer-Verlag, 2002 . |
| [14] | Akaike H. A new look at the statistical model identification. IEEE Trans. Autom. Control , 1974, 19(6): 716-723. DOI:10.1109/TAC.1974.1100705 |
| [15] | Liu L B, Arcone S A. Propagation of radar pulses from a horizontal dipole in variable dielectric ground: A numerical approach. Subsurface Sensing Technologies and Applications , 2005, 6(1): 15-24. |
| [16] | Aki K, Richards P G. Quantitative Seismology. (2nd ed). Sausalito, California: University Science Books, 2002 . |
| [17] | 中国地震局监测预报司. 一九九八年张北地震. 北京: 地震出版社, 1998 : 2 . Monitoring and Forecasting Department of China Earthquake Administration. 1998 Zhangbei Earthquake (in Chinese). Beijing: Seismological Press, 1998 : 2 . |
| [18] | 黄建平, 倪四道, 傅容珊, 等. 综合近震及远震波形反演2006文安地震(Mw5.1)的震源机制解. 地球物理学报 , 2009, 52(1): 120–130. Huang J P, Ni S D, Fu R S, et al. Source mechanism of the 2006 Mw5.1 Wen'an Earthquake determined from a joint inversion of local and teleseismic broadband waveform data. Chinese J. Geophys. (in Chinese) , 2009, 52(1): 120-130. DOI:10.1002/cjg2.v52.1 |
| [19] | Ma S, Custódio S, Archuleta R J, et al. Dynamic modeling of the 2004 Mw6.0 Parkfield, California, earthquake. J. Geophys. Res. , 2008, 113: B02301. DOI:10.1029/2007JB005216 |
| [20] | Mendoza C, Hartzell S. Finite-fault analysis of the 2004 Parkfield, California, earthquake using pnl waveforms. Bull. Seismol. Soc. Am. , 2008, 98(6): 2746-2755. DOI:10.1785/0120080111 |
| [21] | Custódio S, Liu P C, Archuleta R J. The 2004 Mw6.0 Parkfield, California, earthquake: Inversion of near-source ground motion using multiple data sets. Geophys. Res. Lett. , 2005, 32(23): L23312. DOI:10.1029/2005GL024417 |
| [22] | McEvilly T V, Bakun W H, Casaday K B. The Parkfield, California, earthquakes of 1966. Bull. Seismol. Soc. Am. , 1967, 57(6): 1221-1244. |
| [23] | Wu F T. Parkfield earthquake of June 28, 1966: magnitude and source mechanism. Bull. Seismol. Soc. Am. , 1968, 58(2): 689-709. |
| [24] | 梁富康, 于兴河, 李先平, 等. 冀中坳陷深县凹陷的生长断层特点及其对沉积的控制作用. 中国地质 , 2011, 8(2): 263–270. Liang F K, Yu X H, Li X P, et al. Growth faults in Shenxian depression and their control over the sedimentation. Geology in China (in Chinese) , 2011, 8(2): 263-270. |
| [25] | 杨旭升, 刘池阳, 杨斌谊, 等. 冀中坳陷衡水转换断裂带特征及演化. 煤田地质与勘探 , 2004, 32(3): 5–8. Yang X S, Liu C Y, Yang B Y, et al. Characteristics and evolution of Hengshui transfer fault zones in Jizhong depression. Coal Geology & Exploration (in Chinese) , 2004, 32(3): 5-8. |
| [26] | 龙峰, 闻学泽, 徐锡伟. 华北地区地震活断层的震级-破裂长度、破裂面积的经验关系. 地震地质 , 2006, 28(4): 511–535. Long F, Wen X Z, Xu X W. Empirical relationships between magnitude and rupture length, and rupture area, for seismogenic active faults in North China. Seismology and Geology (in Chinese) , 2006, 28(4): 511-535. |
| [27] | 张秀梅. 河北省地震资料汇编. 北京: 地震出版社, 1990 : 324 -327. Zhang X M. Hebei Earthquake Material Collection (in Chinese). Beijing: Seismological Press, 1990 : 324 -327. |