引言

1936年4月1日，广西灵山县东北部(当时属于广东省)发生M6 ¾地震(以下简称“灵山主震”).据文献记载(陈国达，1939)和实地考察结果，当时我国两广各地、香港及越南海防等地均受波及，有感面积达50万km².处于极震区的灵山县和合浦县北5800余间房屋倒塌，200余人受伤，92人死亡.这次地震灾情严重，经济损失惨重(陈国达，1939；广西地震局历史地震小组，1982；顾功叙，1983；谢毓寿，蔡美彪，1985).灵山主震是广西有史以来最大的一次内陆地震，我国研究人员对此作过许多研究(黄河生等，1990；李伟琦，1992；滕吉文等，2001；周本刚等，2008；何军等，2012；任镇寰等，2012)，但是灵山主震的发震时刻、震中等基本参数还需要进一步确定.查看我国不同时期若干版本的地震目录，灵山主震的震中参数出处一致，均为由宏观等震线推算的结果(图 1).大多数地震目录和期刊中，该地震的震中位置为宏观震中(22.5°N，109.4°E)，震级为M6 ¾，而且均未记录发震时刻(李善邦，1960；顾功叙，1983；中国地震局震害防御司，1999；宋志平等，2011；赵荣国等，2012).

图  1  1936年4月1日灵山M6 ¾地震的等震线图

(引自中国地震局震害防御司，1999)

Figure  1.  The isoseismal map of Lingshan M6 ¾ earthquake on April 1, 1936

(after Department of Earthquake Disaster Prevention, China Earthquake Administration, 1999)

下载: 全尺寸图片 幻灯片

据记载，灵山主震发生于1936年4月1日地方时9时31分，相当于格林威治标准时刻2时14分(陈国达，1939；广西地震局历史地震小组，1982).当日2时9分28秒(UTC)，印度尼西亚卡拉克隆岛发生M_W7.7地震，震中为(4.17°N，126.52°E)(Engdahl, Villaseor，2002)，该地震发震时间先于灵山主震不足5分钟，其记录强烈影响或覆盖了灵山地震的记录，造成灵山主震可用数据的稀缺和匮乏，因此，国际上凡是基于仪器记录测定的、应用广泛的地震目录，均未收录灵山主震，例如：20世纪60年代以前，以全球地震观测资料汇编最完整而著称的国际地震汇编(International Seismological Summary, 简写为ISS)，21世纪初Engdahl和Villaseor(2002)汇编的世纪地震目录，Storchak等(2015)近几年完成的国际地震中心全球地震模型(International Seismological Centre-Global Earthquake Model，简写为ISC-GEM)地震目录.基于此，尽管灵山主震已发生80年，但它的基本参数仍然沿用根据地震等震线得出的宏观震中和震级.

本文旨在全球范围内寻找与1936年4月1日灵山主震相关的地震记录，即相应的地震波到时数据和振幅数据，然后根据全球台站的仪器记录，利用现代的参数测定方法和技术对灵山主震进行重新测定，进而提供更具有科学参考意义的基础数据.

1.   资料收集

20世纪30年代，我国仅有少数几个地震台站，还未建立起自己的地震台网.本文所收集的资料时段为1936年4月1日2时至3时(UTC, 下同).收集的资料中，大多数为国外台站的数据，原始资料包括3类：①我国原始的地震记录和文献记载. 1936年，我国仅运行4个台站，即北京鹫峰台，上海徐家汇台，南京水晶台和辽宁营口台，且仅有两个台站的原始记录可供查询，即上海徐家汇台的观测报告、记录图纸和南京水晶台的观测报告，其余两个台站的数据虽被ISS收录，但无法查到原始观测记录；② ISS观测报告中的离散数据. ISS观测报告包括1903—1963年的全球地震观测报告. 1936年4月1日2时至3时，ISS仅编辑了一条地震记录，其发震时刻为1936年4月1日2时9分22秒，震中为(3.6°N，126.7°E)(Richmond，1936)，经查阅该地震的震相到时等数据，证实其与印度尼西亚卡拉克隆岛M_W7.7地震为同一地震. ISS在这条地震后所附的数据为全球各台站的离散数据，即不与任何已知事件相关的数据，本文选取其所记录的全部离散数据，并将各台站的地震波走时数据转换为地震波到时数据；③ ISC-GEM地震目录中与已知地震不相符的数据和离散数据. ISC成立于1964年，其前身为ISS. ISC的主要任务之一为汇编、出版全球最完整的地震观测报告.受GEM资助(ISC, 2012)，ISC组织来自不同国家的学者，自2010年起，历时27个月，完成了ISC-GEM目录，即1900—2009年有仪器记录的全球大震目录(Di Giacomo et al，2015；Storchak et al, 2015), 并在2013年后对该目录进行几次补充和完善. ISC-GEM目录中，1936年4月1日2时至3时也仅有一条地震记录，即上述印度尼西亚卡拉克隆岛M_W7.7地震.本文重点分析了该地震的全球台站记录，将其中明显不属于该地震的震相到时等数据选出，并把这部分数据视为离散数据备用.此外，通过查询ISC-GEM数据库又收录了一些与任何已知事件均无关的离散数据.

综上所述，本文收集到的原始数据集由以上3类数据组成，共收集到全球168个台站的1000余条到时数据，且以第2类数据居多.全球168个台站的分布如图 2所示.

图  2  灵山主震宏观震中和可收集到1936年4月1日2时至3时数据的全球168个台站的分布

Figure  2.  The distrubution of the macro-epicenter of Lingshan mainshock and 168 global stations which could record the information during 02:00--03:00 on April 1, 1936

下载: 全尺寸图片 幻灯片

2.   与灵山主震的数据关联(相关性分析)

不同台站(阵)记录的同一地震的各类震相到时等数据是相关的.将第1节中收集到的全球台站的观测数据与灵山主震进行关联或相关性分析，把属于灵山主震的各台站观测数据组合在一起.

2.1   震源参数初始点范围的选取

为了得到更好的关联结果，在一定范围内不断地进行震源参数初始点的搜索和尝试.以灵山主震的宏观震中(22.5°N，109.4°E)为中心，以100 km为半径所得区域作为初始震中的搜索范围；根据初始震中和台站震中距，查阅走时表获得各台站的初至P波走时，进而得到每个台站所记录到的发震时刻，取多个台站发震时刻的平均值作为初始发震时刻，所取范围为1936年4月1日2时12分至2时14分；初始震源深度范围为0—35 km.

2.2   全球台站的震相数据关联或相关性分析

现代全球台网中，对台站观测数据的关联需考虑4个方面，即初至波到时、S-P波到时差、震级测定和某台站记录到某地震的可能性，用公式表示为

式中：μ为某台站相对于某地震的偏离值，亦为某台站是否属于某地震的定量表示；δt_P为初至波的到时残差，即初至波的观测到时与理论到时之差；δt_S-P为S-P波到时差的残差, 即S-P观测值与理论值之差；δM为震级测定值的误差；l为某台站记录到某地震的可能因素值；ν_S-P, ν_m和ν_l均为常数.

鉴于灵山主震的具体情况和20世纪30年代的观测条件及分析水平，无法对式(1)作出全面分析，仅能考虑其中一项，即初至波的到时残差δt_P，则式(1)可简化为

在计算过程中，预设μ的阈值为μ_max.计算各台站相对于灵山主震的μ值，当μ < μ_max时，该台站的震相等数据即关联到灵山主震中.根据实际情况，不同的台网对μ_max的取值不同，考虑到灵山主震受印尼大地震的干扰，本文将μ_max设定为20 s.

2.3   与灵山主震的关联结果

当最好地满足两个条件时，关联即为成功：①震中距尽可能地连续；②能在关联数据的基础上进行震源参数的修定.经过多次尝试和计算，最终有13个台站与灵山主震关联成功，这些台站的分布如图 3中绿色三角形所示，包括中国南京水晶台、上海徐家汇台、北京鹫峰台、日本丰田台、印度阿格拉台、孟买台、斯洛文尼亚卢布尔雅那台、南非开普敦台、加拿大渥太华台、美国伊萨卡台、福特汉姆台(纽约附近)、乔治城台和橡树岭台.在收集到的全球台站数据中，上述13个台站均有灵山主震的记录.

图  3  灵山主震震中和相应的全球台站分布

Figure  3.  Location of the epicenter of Lingshan mainshock and distribution of the corresponding global stations used for location

Green triangles stand for 13 global stations associated to Lingshan mainshock successfully, and blue triangles for seven global stations used for the relocation with the routine location method and model of China National Seismic Network Center, and red triangles for eight global stations used for the relocation with ISCloc

下载: 全尺寸图片 幻灯片

3.   灵山主震震源参数的测定

灵山主震的数据特点显示，其适合使用大台网的定位方法和一维走时模型.利用上述关联成功的13个全球台站的数据，采用我国国家测震台网常规的定位方法和技术及ISC的定位方法和技术(ISCloc)对灵山主震的震中参数等进行测定.

3.1   采用我国国家测震台网常规的定位方法和技术测定灵山主震

我国国家测震台网常规的定位方法和模型为Geiger法(Geiger, 1912)和J-B走时模型(Bullen, 1963).从20世纪80年代中期至今，中国地震局地球物理研究所出版的《中国地震台临时报告》和《中国地震年报》等报告、中国地震台网中心出版的《中国数字地震台网观测报告》及2009年后产出的统一编目中的M_S≥4.0地震，其日常定位一直沿用Geiger法和J-B走时模型，因而基于该方法和模型产出的地震目录在我国已有30年的系统性和连贯性(陈培善，李强，1989；陈培善等，1993；李保昆等, 2004, 2015).按照我国国家测震台网数据处理的约定，对灵山主震的定位只选用到时残差满足-5 s≤δt_P ≤5 s的初至P波震相，因而上述13个关联成功的台站中仅有7个台站参与定位，即中国南京水晶台、上海徐家汇台、北京鹫峰台、日本丰田台、印度阿格拉台、南非开普敦台和美国橡树岭台, 其分布如图 3中蓝色三角形所示.

使用上述台站数据对灵山主震重新定位后得到其发震时刻为2时12分(20.1±3.64)秒；震中为(22.61°N±45.08 km，109.34°E±32.86 km)；震源深度为(9±27.39 km).震级采用国家台网的测算方法确定为M_S7.0.

3.2   采用ISC的定位方法和技术(ISCloc)测定灵山主震

2011年，ISC发表了其最新的定位方法和技术ISCloc (Bondr, Storchak, 2011)，该方法具有以下特点：使用所有可基于AK135走时模型(Kennett et al, 1995)预测的震相数据；可通过领域算法(neighbourhood algorithm，简写为NA)得到合适的初始点；在使用线性迭代最小二乘法时考虑了走时误差的互相关性；测定自由深度值前需满足一定条件.

用改进的Geiger法(Bondr, McLaughlin, 2009a)求解方程组

式中: G为N×M阶矩阵, 包含N个观测数据对M个变量的偏导数；m为M×1阶震源参数的修正量[ΔT, Δx, Δy, Δz]^T; d为N×1阶走时残差; W为N ×N阶投影矩阵或加权矩阵.

采用奇异值分解(singular value decomposition，简写为SVD)方法求解式(3), 得到

式中：V_W为M×M阶正交矩阵，也称为模型空间；Λ_W为N×M阶对角矩阵；U_W为N×N阶正交矩阵，也称为数据空间.每次迭代中m的改正量为

经过j次迭代，满足精度要求时，有

得到收敛解之后，则定位结果的不确定性(或评价)可由

求得, 式中C_M为模型的协方差矩阵.由式(7)求得一四维椭球体，其投影为误差椭圆、震源深度误差和发震时刻误差.

ISCloc方法中的定位模型为AK135模型，其与J-B走时模型最显著的差别体现在上地幔部分和地核部分.

根据2.3节中的关联结果，使用ISCloc方法可有8个台站参与定位，包括中国上海徐家汇台、北京鹫峰台、日本丰田台、印度阿格拉台、斯洛文尼亚卢布尔雅那台、加拿大渥太华台、南非开普敦台和美国橡树岭台，其分布如图 3中红色三角形所示；除卢布尔雅那台使用S波到时外，其余台站均使用P波到时.对灵山主震重新定位后得到发震时刻为2时12分(17.88±1.88)秒；震中为(22.37°N±47.6 km，109.70°E±22.0 km)；震源深度为默认深度10 km.

4.   讨论与结论

本文收集了1936年4月1日2时至3时全球台站的仪器记录，将其与灵山主震进行关联，分析得到与灵山主震相关联的仪器观测数据.在此基础上，利用现代参数测定方法和技术测定灵山主震发震时刻、震中、震源深度和震级等参数.本文在灵山主震基础参数测定方面取得了新的进展：

1) 利用收集到的全球168个台站1936年4月1日2时至3时的数据，根据宏观震中并通过计算各台站初至P波的到时残差，将其与灵山主震进行关联或相关性分析，得到全球共有13个台站的数据可与灵山主震相关联.

2) 本文利用我国国家测震台网常规的定位方法和技术和ISC的定位方法和技术(ISCloc)对灵山主震进行了参数测定.两种方法和技术的主要差别列于表 1，所测定的灵山主震的结果如表 2和图 4所示.

表  1  CENC和ISCloc的主要差别

Table  1.  The main differences of location methods of CENC and ISCloc

方法	走时模型	定位所用数据类型	走时误差
CENC	J-B	初至P波震相	不相关(各自独立)
ISCloc	AK135	多种震相	相近的传播路径上相关
注：CENC表示我国国家测震台网常规的定位方法和技术，ISCloc表示ISC的定位方法和技术，下同.

下载: 导出CSV 

| 显示表格

表  2  两种方法和技术对灵山主震参数测定的结果

Table  2.  The parameters of Lingshan mainshock determinated by two different methods and techniques

方法	发震时刻 时：分：秒	发震时刻误差/s	走时均方差/s	震中			震中误差椭圆			震源深度 /km	深度误差 /km	台站数	Gap /°	最小震中距 /°	最大震中距 /°	MS
				北纬/°	东经/°		长半轴 /km	短半轴 /km	长轴走向/°
CENC	02:12:20.10	3.64	3.49	22.61	109.34		45.1	32.9	0	9	27.4	7	185	12.60	120.40	7.0
ISCloc	02:12:17.88	1.88	3.64	22.37	109.70		47.6	22.0	162	10.0f	-	8	185	13.66	120.54	-
注：深度数据中的“f”表示默认(缺省)值，Gap为台站最大空隙角.

下载: 导出CSV 

| 显示表格

图  4  灵山主震宏观震中及CENC和ISCloc重新定位后的灵山主震震中

Figure  4.  The macro-epicenter of Lingshan mainshock (green dot), and the epicenters of Lingshan mainshock relocated by CENC (orange dot) and ISCloc (blue dot)

下载: 全尺寸图片 幻灯片

两种方法的定位结果不同，重新定位后的两次震中相距约47.6 km.我国国家测震台网常规的定位方法和技术仅用初至P波到时数据，而且为保证定位精度，P波到时残差大于5 s的台站数据不参与定位计算. ISCloc中，根据IASPEI地面实测数据集中参考地震计算的残差(Bondár, McLaughlin, 2009b; 李保昆等，2014)，震相类型、震中距等因素设置不同的先验测量误差值，并且定位过程中各种震相均有相应的权重.例如，一般地震初至P波、S波和PKP波的先验测量误差值分别设为0.8，1.5，1.3 s，当震相走时残差大于6倍先验测量误差值时，该震相不再参与定位计算.

在灵山主震的深度测定方面，我国国家测震台网常规的定位方法和技术对所用数据无特别约束，得到的是自由深度值；而用ISCloc得到的是默认深度值；在定位过程中固定深度为10 km. ISC对深度测定所用的数据有较严格的规定，当满足以下4个条件中的任何一个时，便计算自由深度值，否则根据AK135模型和震中所在地区给出默认深度值：①至少有一条记录所在台站的震中距在0.2°以内；②至少有来自两个机构的5个以上的深度震相记录，如pP等；③至少有来自两个机构的5个以上的外核面反射震相记录，如PcP, ScS等；④有震中距5°以内的5个以上的S-P震相对.灵山主震不具备上述条件，因而ISCloc不计算其深度，只给出默认深度10 km.

对比两种方法可以看出，ISCloc方法基于现今全球台网观测进行定位，因其可获得丰富的震相，尤其是大量的远震记录，某些地区还有密集的台站分布等；我国国家测震台网常规的定位方法和技术对震相的使用条件相比ISCloc要简单，更适合灵山主震的数据特点，因而首选该方法获取的定位结果.

3) 关于震级的测定.本文收集到的全球168个台站的记录中，均无周期和振幅数据，但上海徐家汇台保存了当时的记录图纸(图 5)，因此本文通过量取徐家汇台的面波周期和振幅，测定震级M_S. 1936年，徐家汇台的水平向观测采用质量为1200 kg的维歇尔地震仪，放大倍率为150，周期为7 s；垂直向观测采用质量为80 kg的维歇尔地震仪，放大倍率为40，周期为6 s；记录方式为熏烟记录.徐家汇台的记录中高频成分覆盖低频成分现象突出，因此量取周期和振幅值较困难，故本文选择三分向记录图中相对较好的东西向记录量取周期T和振幅值A，得到T≈5 s，A≈187 μm.

图  5  1936年4月1日徐家汇台的记录图纸

Figure  5.  The seismograph map of Zikawei on April 1, 1936

下载: 全尺寸图片 幻灯片

ISCloc中，当振幅数据少于3个时，不计算该地震相应的震级.因此本文选用国家测震台网测定震级的方法(李保昆等，2004)来计算灵山主震的震级.

式中Δ为震中距，根据式(8)，(9)和徐家汇台的A, T数据，确定所测震级为M_S7.0.

4) 根据灵山主震的数据特点和20世纪30年代的地震监测水平，无论用什么方法，所测参数均无法达到现代监测条件下的高精度和可信度.从表 2中看出，重新定位后的震中和深度误差均在20 km以上, 震级仅根据徐家汇台的记录所得，单台的结果不足以取代早期的震级M6 ¾.

在数据的收集过程中，本文所用数据与某些台站记录，如越南Phu-lien台和南京水晶台的观测报告等存在矛盾之处，甚至出现一组以上数据关联成功的现象，但是本文首选地震记录图支持的数据组.根据目前的工作条件，无法考证和评价20世纪30年代世界各国的台站数据质量，因而还无法解释某些原始数据之间的矛盾.基于全球台站的仪器记录数据，经上述现代参数测定方法和技术重新测定后得到了新的结果，是灵山主震基础参数测定的重要进展.灵山主震的可用信息掩盖于另一大地震之中，提高对这类历史地震参数的测定精度，将有利于对后续地震活动性、发震构造等研究的深入.