The South Yellow Sea Large Earthquake on June 22, 1518 at Offshore to West of Seoul, South Korea
-
摘要: 1518年6月22日韩国首尔以西海域——南黄海发生大震。地震波及朝鲜半岛全境,并在韩国首尔等沿海地区造成破坏,首尔的烈度为Ⅷ度,余震持续一个多月。地震还影响到中国东部地区。震级定为M7½级。此次地震可能与朝鲜半岛西缘断裂带(南黄海东缘断裂带)的活动有关。震中定于该断裂带附近(36.5°N,125.2°E)。
-
关键词:
- 1518年6月22日地震 /
- 首尔地震 /
- 南黄海大震
Abstract: A large earthquake occurred at offshore to west of Seoul, Korea-South Yellow Sea on June 22, 1518. The shock affected the whole Korean Peninsula which Seoul and other coastal areas and caused damage with the estimated intensity of Ⅷ (MMI) in the Seoul area. Aftershocks last for over a month. The eastern China earthquake was also affected by this earthquake. The magnitude of the earthquake is estimated as M7½.The activities of western margin fault zone of the Korean peninsula (eastern margin fault zone of the Yellow Sea in China) may be responsible for the earthquake. The epicenter of the earthquake is located near the fault zone, at 36.5°N, 125.2°E. -
引言
地震是一种极具破坏性的灾难,全世界地震多发的国家都因地震造成不可弥补的损失,人民生命与财产安全受到严重威胁。对建筑结构采取抗震设防措施是减轻地震灾害损失的根本手段,而地震区划是工程抗震设防的基础和前提条件。现阶段超过设防水平的大地震的发生不能忽视,如唐山地震(1978年)、汶川地震(2008年)等多次地震中,地震实际烈度往往超出抗震设计规范的设防烈度(陈肇元,2008),造成大面积楼房倒塌,数百万人流离失所。因此,灾前对既有建筑结构开展震害评估,并采取相应防灾措施尤为重要。传统的地震预测分析方法通常指对结构在特定地震作用下发生不同破坏状态的可能性进行分析,由于未充分考虑地震活动的不确定性,不能较好地反映建筑结构潜在的地震风险水平。Preciado等(2015)进行建筑结构震害风险评估时联合结构易损性和场地危险性,并以二者卷积为表征,表示结构在未来某一段时间内超过或发生某种地震破坏状态的概率,从而明确度量建筑结构未来可能发生的风险水平,为做好灾前预防、地震应急准备和灾后救援工作提供科学的决策依据,也是贯彻“从减少灾害损失向减轻灾害风险转变”(罗华春,2017)思想的体现。
在长期的实践探索中,国内外学者在结构震害风险研究方面取得了一系列有价值的成果。Melani等(2016)采用IDA方法对某低层钢筋混凝土框架结构开展易损性分析,并基于框架结构概率风险四次积分公式,对其进行地震风险分析和年损失估计。Berto等(2012)基于强震环境对钢筋混凝土构件的影响研究,提出概率风险管理方法,计算与建筑结构长期性能相关的地震风险。谷音等(2012)针对某典型矮塔斜拉桥,采用LHS-MC方法对其地震风险概率进行评估。尹超等(2016)以实际路堤结构为例,采用蒙特卡罗抽样方法评价有无挡土墙路堤的震害风险,提出基于风险可接受度的路堤抗震设计与管理方法。然而,城市中存在的大量钢筋混凝土框架结构震害风险尚未可知,对不同抗震设防烈度区框架风险等级的划分仍处于研究阶段。
基于上述认识,综合考虑我国地域辽阔且国土范围内地震危险性差异大的实际情况(马玉宏等,2002)和目前国内最主要的钢筋混凝土建筑结构形式,本文以钢筋混凝土框架结构为例,按照不同设防烈度要求设计3个结构模型,采用IDA方法对模型进行地震易损性分析和不同地震作用的可能性分析,评估结构震害风险,为抗震防灾对策的制定提供参考。
1. 结构易损性分析
根据地震影响区域、地震损伤数据采集方法和计算方法的不同,地震易损性评价的方法主要有经验方法和理论方法(Wu等,2012)。理论方法内容主要包括确定损伤判别条件、确定结构模型、选择和输入地震动、IDA分析、绘制易损性曲线等,主要通过计算数值反映结构各损伤状态与地震动强度之间的概率关系(孙柏涛等,2018)。IDA方法是目前性能化地震工程(Performance-Based Earthquake Engineering,PBEE)中评价结构地震反应、预测结构震害最常用的方法之一。该方法可对结构进行从弹性状态到弹塑性状态直至倒塌的全过程分析,能较好地反映结构在未来不同地震作用下的抗震能力变化规律(侯炜,2013)。因此,本文采用IDA方法进行结构地震易损性分析。
1.1 IDA指标选取
IDA分析结果由地震动强度指标IM(Intensity Measures)与结构性能参数指标DM(Damage Measures)之间的关系表示,目前国内外常用的地震动强度指标主要包括第一模态谱加速度Sa(T1,5%)和地震加速度峰值PGA(g),为了与地震危险性关联,本文选取PGA作为地震动强度指标。最大层间位移角${\theta _{\max }}$通过体现层间变形的综合结果在整体层面反映结构抗震性能(蒋欢军等,2009),因此,本文选取${\theta _{\max }}$作为结构性能参数指标。
1.2 结构破坏等级判别
建筑结构地震破坏等级随着破坏程度的加深分为基本完好、轻微破坏、中等破坏、严重破坏和倒塌。框架结构最大层间位移角${\theta _{\max }}$与地震破坏等级之间的对应关系如表 1所示(马东辉等,2007)。
表 1 最大层间位移角与地震破坏等级的关系Table 1. The relationship between the maximum inter-story drift ratios and damage degrees破坏等级 基本完好 轻微破坏 中等破坏 严重破坏 倒塌 ${\theta _{\max }}$ ${\theta _{\max }}$<1/400 1/400≤${\theta _{\max }}$<1/250 1/250≤${\theta _{\max }}$<1/125 1/125≤${\theta _{\max }}$<1/50 ${\theta _{\max }}$≥1/50 1.3 结构模型信息
本文以住宅中常用的内廊式结构为例进行研究,考虑所建模型应能从一定程度上反映按照不同设防标准设计和建造的一般框架结构抗震性能,采用较规则的平面布置,如图 1所示。按照相关规范(中华人民共和国住房和城乡建设部等,2010;中华人民共和国住房和城乡建设部,2012)要求,采用PKPM软件建立3个抗震设防烈度分别为6度(0.05g)、7度(0.1g)、8度(0.2g)的3跨10层框架结构模型。场地类别为Ⅱ类,设计地震分组为第二组,场地特征周期为0.40s。结构首层层高4.5m,其他层层高3.6m,建筑高度和抗震等级如表 2所示。恒荷载标准值取4.6kN/m2,活荷载标准值不上人屋面取0.5kN/m2、楼面取2.5kN/m2。
表 2 结构抗震等级Table 2. Seismic grade of structures项目 模型1 模型2 模型3 设防烈度 6度(0.05g) 7度(0.1g) 8度(0.2g) 建筑高度/m 36.9 36.9 36.9 抗震等级 三级 二级 一级 1.4 地震动选取和输入
不同地震动输入得到的地震作用效应值差异较大,有时可达数十倍(王亚勇等,1991)。如果建筑物所在场地具有实际强震记录,则为最优选择,但大部分情况下并不具备满足条件的强震记录。与美国地质勘测中心(United States Geological Survey,USGS)划分的S2场地基本相似,为使选用的地震动频谱特性尽可能与建筑物所在场地保持一致,从S2场地中选取30条具有代表性的强震记录(FEMA建议20条以上),将PGA调幅后作为地震动输入,如表 3所示。
表 3 地震动基本信息Table 3. Basic information of ground motion编号 名称 时间/年 台站 PGA/g 1 Livermore-02 1980 San Ramon-Eastman Kodak 0.191 2 Westmorland 1981 Parachute Test Site 0.219 3 Imperial Valley-06 1979 Delta 0.284 4 Northridge-06 1994 LA-Century City CC North 0.123 5 Whittier Narrows-01 1987 Downey-Co Maint Bldg 0.177 6 Whittier Narrows-01 1987 Brea Dam(Downstream) 0.231 7 Whittier Narrows-02 1987 LA-116th St School 0.157 8 Imperial Valley-02 1940 El Centro Array #9 0.342 9 Imperial Valley-06 1979 El Centro Array #7 0.42 10 Northridge-01 1994 LA-Century City CC North 0.223 11 Northridge-01 1994 LA-Obregon Park 0.467 12 Sierra Madre 1991 LA-Obregon Park 0.224 13 Superstition Hills-02 1987 Parachute Test Site 0.451 14 Loma Prieta 1989 Capitola 0.480 15 Victoria,Mexico 1980 LA-Obregon Par 0.118 16 Morgan Hill 1984 Capitola 0.117 17 San Fernando 1971 Carbon Canyon Dam 0.070 18 San Fernando 1971 Gormon-Oso Pump Plant 0.087 19 Imperial Valley-07 1979 Bonds Corner 0.092 20 Imperial Valley-07 1979 El Centro Array #7 0.161 21 Livermore-01 1980 San Ramon-Eastman Kodak 0.106 22 Victoria,Mexico 1980 Chihuahua 0.118 23 Whittier Narrows-02 1987 Inglewood-Union Oil 0.133 24 Landers 1992 Downey-Co Maint Bldg 0.046 25 Landers 1992 Inglewood-Union Oil 0.040 26 Northridge-06 1994 LA-Baldwin Hills 0.052 27 Hector Mine 1999 Downey-Co Maint Bldg 0.031 28 Hector Mine 1999 LA-116th St School 0.022 29 CA/Baja Border Area 2002 Calexico Fire Station 0.094 30 CA/Baja Border Area 2002 El Centro Array #7 0.078 1.5 易损性分析
统计和实际工程中认为结构能力参数C服从对数正态分布,假定结构反应参数D也服从对数正态分布(李玉,2015):
$$C = \ln (\mathop C\limits^ \wedge, {\beta _{\rm{c}}})$$ (1) $$D = \ln (\mathop D\limits^ \wedge, {\beta _{\rm{d}}})$$ (2) 式中,$\mathop C\limits^ \wedge $、$\mathop D\limits^ \wedge $分别表示结构能力和结构反应均值;${\beta _{\rm{c}}}$、${\beta _{\rm{d}}}$分别表示结构能力和结构反应对数标准差。
结构进行IDA分析后,可得最大层间位移角${\theta _{\max }}$(结构反应D)和$PGA$的对应关系(何益斌等,2016)为:
$$D = aPG{A^b}$$ (3) 上式两边分别取对数可得:
$$\ln D = \ln a + b\ln PGA = A + B\ln PGA$$ (4) 式中,a、b、A、B为回归系数,可根据结构IDA分析结果拟合得到。
易损性曲线表示不同地震作用下结构反应D超出某一破坏等级所定义的结构能力C的条件概率(贾晗曦等,2019),即结构失效概率${P_{\rm{f}}}$:
$${P_{\rm{f}}} = P(C/D \leqslant 1)$$ (5) 根据式(1)和式(2),已知D和C均服从对数正态分布,则结构失效概率可表示为:
$${P_{\rm{f}}} = \mathit{\Phi} \left[ {\frac{{\ln (\mathop D\limits^ \wedge /\mathop C\limits^ \wedge)}}{{\sqrt {\beta _{\rm{c}}^{\rm{2}} + \beta _{\rm{d}}^2} }}} \right]$$ (6) 式中,$\mathop C\limits^ \wedge $根据表 1中以最大层间位移角表示的结构性能量化指标取值;$\sqrt {\beta _d^{\rm{2}} + \beta _c^{\rm{2}}} $根据规定:当结构易损性曲线以地震峰值加速度$PGA$为自变量时取0.5,本文取0.5;$\mathit{\Phi}( \cdot )$为正态分布函数,其取值可通过查询标准分布表确定。
对30条地震动PGA进行调幅,分别调至0.2g、0.4g、0.6g、0.8g、1.0g、1.2g,得到结构性能参数(θmax)与地震动强度(PGA)之间的关系曲线,即IDA曲线簇,如图 2—4所示。
对结构模型进行IDA分析后,将横纵坐标取对数,拟合线性关系,并根据上述计算方法分别得到结构易损性曲线,如图 5—7所示。
由图 5—7可知,结构发生各种破坏状态的累积超越概率曲线均随着PGA的增大呈上升趋势,PGA越大,结构震害越严重;对于同一破坏状态的累积超越概率曲线,8度设防烈度结构上升趋势最缓,同一地震动强度作用下8度设防烈度结构发生严重破坏和倒塌的概率更小,说明8度设防烈度结构较6度和7度设防烈度结构具有更好的抗震能力。
2. 震害风险评价
2.1 地震烈度概率分布
根据高小旺等(1985)对华北、西北、西南地区的45个城镇地震危险性分析结果,一个地区在未来50年内发生的地震烈度服从极值Ⅲ型分布,分布函数为:
$$F(x) = \exp \left[ { - {{\left({\frac{{\omega - x}}{{\omega - \varepsilon }}} \right)}^K}} \right]$$ (7) 对分布函数求导可得概率密度函数$f(x)$,即某一场地未来50年内发生各地震烈度的概率密度:
$$f(x) = \frac{{K{{(\omega - x)}^{K - 1}}}}{{{{(\omega - \varepsilon)}^K}}}\exp \left[ { - {{\left({\frac{{\omega - x}}{{\omega - \varepsilon }}} \right)}^K}} \right]$$ (8) 式中,$x$表示地震烈度,为1—12的离散变量;$\omega $为烈度上限值,根据目前通用的地震烈度划分方法,$\omega $取12;$\varepsilon $为众值烈度,即未来50年内超越概率为63.2%的地震烈度;K为形状参数,取决于一个地区地震背景的复杂性,根据文献(吕大刚,1999),形状参数K的取值如表 4所示。由式(8)绘制不同设防烈度地区地震概率密度曲线,如图 8所示。
表 4 形状参数K的取值Table 4. Values of shape parameter K基本烈度 6度 7度 8度 K 9.7932 8.3339 6.8713 2.2 震害风险计算方法
结构震害风险可定义为:在考虑场地地震危险性的基础上,结构发生各种破坏状态的可能性,在数值上等于场地危险性和地震易损性的卷积。结构在未来50年内发生震害的概率P计算公式如下:
$$P = \sum\limits_{{\rm{all}}\;{x_i}} {P(C \leqslant D|PGA = {x_i})} P(PGA = {x_i})$$ (9) 式中,${x_i}$为未来可能遭遇的$PGA$。由于$PGA$的取值是连续的,因此式(9)可写为:
$$P{\rm{ = }}\int {P(C \leqslant D|PGA = {x_i})} f({x_i}){\rm{d}}x$$ (10) 式中,$f({x_i})$为$PGA$取值概率密度函数,与地震烈度概率密度函数$f(x)$一致。
2.3 基于蒙特卡罗方法的风险评价
蒙特卡罗方法是以概率和统计理论方法为基础的数值计算方法,也称统计模拟方法或随机抽样技术,通过使用随机数(或伪随机数)解决计算问题。本文根据蒙特卡罗思想产生大量符合地震概率密度的随机数,模拟未来场地可能遭遇的地震烈度。考虑积分计算的精确性,将地震烈度作为连续变量,以式(7)为母体,采用蒙特卡罗方法通过MATLAB软件对地震烈度进行n次抽样,并根据公式将地震烈度转化为PGA,则每个PGA在整体中出现的概率为$1/n$,产生的随机样本n越大,结果越接近真实值。本文取n=50000,50000个随机数即表示该地区可能发生的地震烈度值及其分布特征。当对不同设防烈度地区进行计算时,由于其分布特征不同,因而对结构造成的危险性也不同。
将抽样后得到的地震烈度I根据刘恢先等(1981)提出的公式转换为PGA:
$$ PGA = {10^{\left( {I \cdot \lg 2 - 0.01} \right)}} $$ (11) 根据大数定律,式(10)可写为:
$$P = \frac{1}{n}\sum\limits_{i = 1}^{50000} {P(C \leqslant D|PGA = {x_i})} $$ (12) 结构超越各级破坏的风险概率${P_j}$为:
$${P_j} = \sum\limits_{i = 1}^{50000} {\frac{{P(C \leqslant D|PGA = {x_i})}}{{50000}}} $$ (13) 式中,j=1、2、3、4、5,代表每级震害状态,分别对应基本完好、轻微破坏、中等破坏、严重破坏、倒塌。
同理,可知结构发生各级震害的风险${C_j}$为:
$${C_j} = \left\{ \begin{gathered} {P_j} - {P_{j + 1}}\;\;\;\;\;\;\;j \leqslant 4 \\ {P_j}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;j = 5 \\ \end{gathered} \right.$$ (14) 2.4 震害风险评价
依据上述分析模型,结合本文结构易损性分析结果,计算得到不同设防烈度框架结构50年震害风险结果,如表 5—7所示。
表 5 6度设防烈度框架结构震害风险评价结果(%)Table 5. Seismic risk assessment results of 6 intensity frame structure(%)震害等级 基本完好 轻微破坏 中等破坏 严重破坏 倒塌 超越风险 100 29.04 8.91 1.92 0.14 发生风险 70.96 20.13 6.99 1.78 0.14 表 6 7度设防烈度框架结构震害风险评价结果(%)Table 6. Seismic risk assessment results of 7 intensity frame structure(%)震害等级 基本完好 轻微破坏 中等破坏 严重破坏 倒塌 超越风险 100 42.34 16.09 4.25 0.43 发生风险 57.66 26.25 11.84 3.82 0.43 表 7 8度设防烈度框架震害风险评价结果(%)Table 7. Seismic risk assessment results of 8 intensity frame structure(%)震害等级 基本完好 轻微破坏 中等破坏 严重破坏 倒塌 超越风险 100 55.20 25.64 8.28 1.07 发生风险 44.80 29.56 17.36 7.21 1.07 由表 5—7可知,在3类设防烈度地区结构处于基本完好和轻微破坏状态风险的概率较大,分别为91.09%、83.91%和74.36%,发生严重破坏和倒塌状态风险的概率较低。由6度、7度、8度抗震设防烈度框架结构震害风险对比可知,8度设防烈度框架结构发生严重破坏和倒塌状态风险的概率最大,分别为7.21%、1.07%;7度设防烈度框架结构次之,发生严重破坏和倒塌状态风险的概率分别为3.28%、0.43%;而6度设防烈度框架结构发生严重破坏和倒塌状态风险的概率最低,分别为1.78%、0.14%。综合结构地震易损性结果可以看出,造成结构震害风险的因素除结构抗震能力外,还与其所处抗震设防烈度地区有关。
2.5 震害风险等级划分
为综合考虑结构不同震害等级的发生风险,并实现对不同设防烈度地区结构震害风险的排序,本文采用综合震害风险指数R表征结构震害风险平均程度,其数值范围0—1,结构综合震害风险指数越小,表示结构震害风险水平越低。参考震害指数(尹之潜等,1990)划分方法,建议各震害等级对应的综合震害风险指数范围如表 8所示。
表 8 综合震害风险指数划分Table 8. Composite seismic risk indexes classification震害等级 基本完好 轻微破坏 中等破坏 严重破坏 倒塌 综合震害风险指数范围 0<R1≤0.2 0.2<R2≤0.4 0.4<R3≤0.6 0.6<R4≤0.8 0.8<R5≤1 综合震害风险指数中值 0.1 0.3 0.5 0.7 0.9 结构综合震害风险指数计算公式如下:
$$R = \sum\limits_{i{\rm{ = 1}}}^5 {\sum\limits_{j = 1}^5 {{R_i}{C_j}} } $$ (15) 式中,Cj为结构发生各级震害的风险;Ri为各级震害风险指数中值,i=1、2、3、4、5。
结构综合震害风险指数由式(15)计算,既考虑结构各级震害风险的发生概率,又量化结构震害风险。
为了更直观地对建筑结构震害风险水平进行评定和分类,本文根据结构综合震害风险指数将结构震害风险等级划分为4级,由低到高分别为Ⅰ级、Ⅱ级、Ⅲ级、Ⅳ级,从Ⅰ级到Ⅳ级表示结构震害风险越来越严重,如表 9所示。
表 9 震害风险等级划分Table 9. Seismic risk grades classification综合震害风险指数 0<R≤0.25 0.25<R≤0.5 0.5<R≤0.75 0.75<R≤1 震害风险等级 Ⅰ级 Ⅱ级 Ⅲ级 Ⅳ级 结合表 5—7中不同设防烈度框架结构发生各震害等级的风险,按照式(15)进行计算,得到结构综合震害风险指数,进而根据表 9中综合震害风险指数范围判断结构震害风险等级,计算和判断结果如表 10所示。
表 10 结构综合震害风险指数和震害风险等级Table 10. Composite seismic risk indexes and risk grades of buildings结构 6度设防烈度框架结构 7度设防烈度框架结构 8度设防烈度框架结构 综合震害风险指数 0.18 0.22 0.28 震害风险等级 Ⅰ级 Ⅰ级 Ⅱ级 经过计算,6度、7度、8度设防烈度框架结构综合震害风险指数分别为0.18、0.22、0.28,可知8度设防烈度地区框架结构震害风险最高,7度设防烈度地区框架结构震害风险次之,6度设防烈度地区框架结构震害风险最低。根据表 9可知,6度、7度设防烈度框架结构综合震害风险指数为0—0.25,震害风险等级为Ⅰ级;8度设防烈度框架结构震害风险指数为0.25—0.5,震害风险等级为Ⅱ级。
3. 结论
1)进行建筑结构震害风险分析时,综合考虑结构地震易损性和地震危险性,通过对风险水平进行量化与评估,更全面和直观地反映建筑结构未来一段时间存在潜在破坏可能,对制定科学合理的抗震防灾对策更有意义。
2)在抗震能力方面,按照6度、7度、8度设防烈度设计的框架结构随着设防烈度的升高,抗震性能有所提高。但由于不同设防烈度地区地震危险性的差异,导致综合考虑地震危险性后结构震害风险趋势发生变化,如8度设防烈度地区结构震害风险等级为Ⅱ级,而6度和7度设防烈度地区结构震害风险等级为Ⅰ级。因此地震危险性在震害风险中起着不容忽视的作用,在城市抗震防灾规划中应注意高烈度地区建筑结构震害风险的防范。
3)由于建筑结构使用功能不同,使其平立面布置等具有较大差异,必然导致其在结构设计方面存在差异,难以通过单一模型表示。本文仅以3个3跨10层内廊式规则框架结构为例,评估其在不同设防烈度地区的震害风险,存在一定局限性。建筑结构模型的不确定性、地震活动的不确定性、认知不确定性及其他因素对风险的影响有待进一步研究。
-
顾功叙, 1983.中国地震目录(公元前 1831 年-公元 1969 年).北京:科学出版社. 国家地震局震害防御司, 1995.中国历史强震目录(公元前 23 世纪-公元 1911 年).北京:地震出版社, 365-367 韩国国史编纂委员会, 2016.国史数据库. http://sillok.history.go.kr. 郝天珧, Mancheol S.,王谦身等, 2002.根据重力数据研究黄海周边断裂带在海区的延伸.地球物理学报, 45(3):385-397. http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX200203009.htm 郝天珧,刘建华, Mancheol S.等, 2003.黄海及其邻区深部结构特点与地质演化.地球物理学报, 46(6):803-808. http://www.cnki.com.cn/Article/CJFDTOTAL-DQWX200306012.htm 李钦祖,于利民,姚振兴等, 1994.黄海地震带是一条重要地震带.见:陈运泰主编,中国固体地球理学进展.北京:海洋出版社, 339-345. 李钦祖,于利民,刁桂苓等, 1997.中国地震科学的特色.见:陈运泰主编,中国地震学研究进展——庆贺谢毓寿教授八十寿辰.北京:地震出版社. 李裕澈, 李德基, 吴锡薰等, 2003. 1996 年 11 月 9 日长江口以东海域 MS 6.1 地震对韩国的影响及烈度分布. 地震学报, 25(4):446-448. http://www.cnki.com.cn/Article/CJFDTOTAL-DZXB200304012.htm 刘昌森,景天永,孙庆煊等, 2002.苏·浙·皖·沪地震目录(公元 225-2000 年).北京:地震出版社, 226, 126. 谭其骧, 1982.中国历史地图集:第七册(元·明时期).北京:中国地图出版社. 汪素云,时振梁, 1992.利用地震有感范围判定震级.见:国家地震局地球物理研究所编,中国地震考察(公元前 466 年-公元 1900 年).北京:地震出版社, 104-108. 吴戈,房贺岩,李志田等, 1992.东北地震史料辑览.北京:地震出版社, 17. 吴戈,刘昌森,翟文杰等, 2001.黄海及其沿岸历史地震编目与研究.北京:地震出版社, 34. 谢毓寿,蔡美彪, 1985.中国地震历史资料汇编:第二卷.北京:科学出版社, 251. 鄢家全,张志中,王健等, 2011.中国历史地震烈度表研究.地震学报, 33(4):515-531. http://www.cnki.com.cn/Article/CJFDTOTAL-DZXB201104012.htm -