引言

强震观测在核电厂、高铁等重大工程抗震设防中具有不可替代的作用。全球几乎所有的核电厂都布设有强震观测系统(李山有等,2004)。在核电厂设计时,均应按照相关准则、规范,在特定的位置和设备上设置不同的地震监测仪器,主要用于地震期间为操作员提供信息,从而使其决定应采取的即时应急操作和核电厂进一步的运行方式。2007年7月16日,日本新潟发生6.8级地震,此次地震虽然对核电厂造成强烈影响,但地震触发了安全装置,使4台正在运行的核反应堆自动停堆(陈志高等,2015);2011年3月11日,日本东北部海域发生9.0级强烈地震,共有11台核电机组自动停堆,停堆后反应堆的热功率迅速下降到4%,但由于这4%的余热处理不当,最终引起福岛核泄露事故(何召壮,2012陈志高等,2015)。这2次日本核事故证明,地震监测系统能在强地震中发挥重要作用(陈志高等,2015)。中国的相关核安全法规导则并未提及地震自动停堆的要求,但在《“十二五”期间新建核电厂安全要求(报批稿)》的“要求25:仪控系统”中,增加了“宜设置地震自动停堆系统”(郑华等,2015)。同样,中国正在研究的高铁地震仪表系统,可在发生高烈度地震时,使高铁紧急制动,尽量避免人员伤亡和财产损失。虽然地震仪表紧急停堆、紧急制动系统有助于提高核电厂、高铁的安全性,但如果不能准确地通过地震仪表监测系统评估地震对核电厂、高铁的破坏作用,则有可能发生误报警、误停堆和误停车的现象,将会造成大量不必要的资源浪费,甚至事故。

应用地震动参数快速计算仪器烈度,对于震后快速评估地震的破坏程度、采取相应紧急措施具有重要作用,研究人员对此进行了大量的研究。美国和日本的仪器烈度计算方法已经应用到实际的工程当中,并取得了非常好的效果,减少了财产损失和人员伤亡。美国ShakeMap系统采用峰值加速度(PGA)和峰值速度(PGV)进行仪器烈度计算,最初该系统只处理发生在美国南加州地区的地震,如今则可对全球范围内任何地区发生的破坏性地震快速产出一系列结果(Wald等,1999Brackman等,2006),包括峰值地面运动加速度等值图、峰值地面运动速度等值图、仪器烈度分布图以及反应谱等值图等(金星等,2013)。日本气象厅(JMA)则采用地震动时程三分量、持时、频谱等多个地震动要素进行计算,发布的产品也非常丰富,包括PGA等值图、PGV等值图、气象厅测定烈度分布图、周期1—2s(阻尼比为5%)的平均反应谱等值图、谱烈度等值图以及质点运动轨迹图等(金星等,2013)。

本文将对仪器烈度计算方法进行分析,以期为中国核电厂、高铁地震仪表监测系统的研究提供参考。

1 地震动参数的选取
1.1 各国仪器烈度地震参数

美国ShakeMap系统的仪器烈度综合考虑了幅值和频谱2个参数,在烈度较低的情况下采用峰值加速度(PGA)进行计算,在烈度较高的情况下采用峰值速度(PGV)进行计算,当烈度介于Ⅴ度和Ⅶ度之间时,采用峰值加速度(PGA)和峰值速度(PGV)共同计算仪器烈度。不同的烈度等级采用不同的幅值进行计算,主要考虑了频谱的影响,仪器烈度计算公式(Wald等,1999, 2006)如下:

$ {{I}_{\text{mm}}}=\left\{ \begin{align} &\text{2}.20\times \text{lg}(\text{PGA})+1.00, \ \ \ \ \text{I}\le \text{ }{{I}_{\text{mm}}}<\text{IV} \\ &\text{3}.6\text{6}\times \text{lg}(\text{PGA})-1.66, \ \ \ \ \text{V}\le {{I}_{\text{mm}}}<\text{VIII} \\ &\text{3}.4\text{7}\times \text{lg}(\text{PGV})+2.35, \ \ \ \ \text{V}\le {{I}_{\text{mm}}}<\text{IX} \\ &\text{2}.1\text{0}\times \text{lg}(\text{PGV})+3.40, \ \ \ \ \text{I}\le \text{ }{{I}_{\text{mm}}}<\text{IV} \\ \end{align} \right. $ (1)

相比美国算法,日本气象厅(JMA)的算法增加了持时的影响,在考虑幅值和频谱时也有很大不同。JMA采用滤波的方法考虑频谱对烈度的影响,滤波频带主要根据其建筑特点确定;同时,采用全振动向量的方法考虑幅值对烈度的影响,选用持时$ \ge $0.3s的全振动向量计算仪器烈度值。计算步骤(Kuwata等,2002)为:首先,将3个时程分量$ g(t) $分别作傅里叶变换得到$ {{G}_{i}}(\omega) $;然后,将$ {{G}_{i}}(\omega) $进行滤波处理得到$ {{{G}'}_{i}}(\omega) $,见式(2):

$ {{{G}'}_{i}}(\omega)={{G}_{i}}(\omega)\times {{F}_{1}}(\omega)\times {{F}_{2}}(\omega)\times {{F}_{3}}(\omega) $ (2)

式中,F1为调幅滤波器,F2为低通滤波器,F3为高通滤波器,分别为:

$ {{F}_{1}}(\omega)={{({1}/{\omega }\;)}^{{1}/{2}\;}} $ (3)
$ {{F}_{2}}(\omega)={{(1+0.694{{x}^{2}}+0.24{{x}^{4}}+0.0557{{x}^{6}}+0.009664{{x}^{8}}+0.00134{{x}^{10}}+0.000155{{x}^{12}})}^{{-1}/{2}\;}} $ (4)
$ {{F}_{3}}(\omega)={{(1-\text{exp}{{(-{\omega }/{{{\omega }_{0}}}\;)}^{3}})}^{{1}/{2}\;}} $ (5)

其中,$ x={\omega }/{{{\omega }_{c}}}\; $$ {{\omega }_{c}}=10\text{Hz}$$ {{\omega }_{0}}=0.5\text{Hz}$$ \omega $为地震波的频率。

$ {{{G}'}_{i}}(\omega) $傅里叶逆变换,得到$ {{{g}'}_{i}}(t) $,再将三分向$ {{{g}'}_{i}}(t) $合成向量加速度:

$ A=\sqrt{{{{{g}'}}_{1}}{{(t)}^{2}}+{{{{g}'}}_{2}}{{(t)}^{2}}+{{{{g}'}}_{3}}{{(t)}^{2}}} $ (6)

最后,将全振动时程上持时$\ge $0.3s的加速度($ \tau ({{A}_{m}})\ge 0.3\text{s} $)带入式(7),求得仪器烈度:

$ {{I}_{\text{JMA}}}=2\lg {{A}_{m}}+0.94 $ (7)

袁一凡(1998)曾提出了1套基于模糊方法计算地震烈度的算法,综合考虑频谱、幅值、持时3个要素,选取了8个地震动参数进行计算。

1.2 地震数据的选取

强震数据的选取对于分析结果具有重要影响。为使选取的强震数据尽可能合理,分别从软土、硬土、基岩3种不同场地,收集了1982—2008年间部分国内的强震记录,共计612条,震级范围3.2—8.0级,震中距0—1600km。由于目前对近场没有明确的定义(周正华等,2002),故本文选择不同的震中距范围,分别为0—10km、10—30km、30—100km、100—500km和500km以上。

1.3 滤波频带的选取

日本的算法使用的滤波频带为0.5—10Hz,这主要考虑了日本的房屋建筑的振动特性(袁一凡,1998)。本文滤波频带的选取主要考虑3个方面(李亮,2011):①滤波频带对加速度峰值的影响。统计发现当高频频率截止到10Hz时,加速度的峰值比较稳定,当高于10Hz时,滤波的效果就不明显了,因此建议将高频截止频率定为10Hz。②大震远震、近震小震和不同烈度区加速度记录的频率成份的分布。经统计分析可以发现近震小震高频成份较丰富,远震大震低频成份较丰富。当低频截止频率为0.1Hz时,完全可以包含强震记录的主要频带范围。③考虑中国建筑结构的基本自振频率。根据《建筑结构荷载规范》中的规定,高层建筑结构的基本自振频率经验公式为:钢结构T1=[0.10,0.15]H,钢筋混凝土结构T1=[0.05,0.10]H,假设H取值为50层,那么高层结构的自振频率为0.13—0.4Hz,完全包括在本文所选的0.1—10Hz范围之内。所以本文所选的滤波频带范围比较符合中国建筑的实际情况。通过分析不同频带对加速度时程变化的影响,发现滤波对绝大多数的加速度时程影响不大,变化率小于20%。但对震中距较小的小震,滤波的效果比较明显,有的滤波前后变化率甚至达到了90%,根据统计分析结果并结合中国建筑结构的实际特点,认为0.1—10Hz滤波频带比较合适。

1.4 加速度的选取

袁一凡(1998)共选用了8个地震动参数:包括2个水平向峰值加速度、垂直分向加速度、卓越频率、20%相对持时和8Hz、5Hz、2Hz、1Hz周期点的反应谱值(刘锡荟等,1983冯德益等,1989苏经宇等,1991)。日本采用合成加速度持时≥0.3s的等效峰值加速度,能较好地反映地震对结构的破坏情况。但是选用多大的持时更为合适,仍是值得探讨的问题。李亮(2011)经统计发现,三分向加速度时程合成的全振动时程(简称“全加速度Aall”)和有效峰值加速度有很好的对数线性关系(分别选取了持时≥0.3s的A0.3、持时≥0.5s的A0.5)。选择不同持时的有效峰值加速度可能对结构的破坏具有不同的作用,但是带入对数烈度经验公式后,所得到的仪器烈度值差别不大,因此本文建议直接选择全加速度Aall计算烈度值。

2 推荐的仪器烈度算法

综上所述,本文拟采用0.1—10Hz滤波频带进行滤波,并选取三分向合成全振动(速度)时程带入拟合公式计算仪器烈度值。采用分档计算的方法计算仪器烈度,具体步骤如下:

(1)对三分向地震动时程分别滤波(滤波频带为0.1—10Hz)。

(2)将三分向加速度时程a、速度时程v分别合成加速度时程Aallt)、速度时程Vallt)。

$ {{A}_{\text{all}}}(t)=\sqrt{a_{\text{EW}}^{\text{2}}+a_{\text{NS}}^{\text{2}}+a_{^{\text{UD}}}^{\text{2}}} $ (8)
$ {{V}_{\text{all}}}(t)=\sqrt{v_{\text{EW}}^{\text{2}}+v_{\text{NS}}^{\text{2}}+v_{\text{UD}}^{\text{2}}} $ (9)

(3)将AallVall带入烈度经验公式,计算求得仪器烈度值。

烈度在Ⅸ度以上时仪器烈度值采用全速度拟合的公式计算,Ⅴ度以下时仪器烈度值采用全加速度拟合的公式计算,Ⅵ—Ⅷ度采用全加速度、全速度共同拟合的公式计算。由于有烈度调查的强震台记录主要分布在Ⅵ—Ⅷ度间,故本文只给出了Ⅵ—Ⅷ度拟合的烈度公式:

$ I=\left\{ \begin{align} &\text{1}\text{.205}\times \text{lg}(A_\text{all})+\text{4}\text{.238}, \ \ \ \text{VI}\le I\le \text{VIII} \\ &\text{0}\text{.985}\times \text{lg}(V_\text{all})+\text{5}\text{.87}, \ \ \ \ \ \text{VI}\le I\le \text{VIII} \\ \end{align} \right. $ (10)

拟合的强震数据主要来自5次西部地区的地震,分别为2007年6月3日宁洱6.4级地震、2008年8月31日攀枝花6.1级地震、2009年7月9日姚安6.0级地震、2008年8月21日盈江5.9级地震和2008年10月5日新疆6.8级地震,见表 1李敏,2010)。

表 1 强震记录 Table 1 Strong motion records

本文推荐的仪器烈度计算方法流程为:将三分向时程合成为加速度时程Aall和速度时程Vall;将Vall带入Ⅸ度以上经验公式,若计算得到的烈度值大于等于Ⅸ度,则输出;若计算得到的烈度值介于Ⅵ—Ⅷ度间,则用公式(10)计算,否则,带入Aall经验公式,具体流程如图 1所示。需要说明的是,正如前文所述,由于本文所收集到的有烈度调查的强震台记录主要分布在Ⅵ—Ⅷ度之间,因此这里并没有拟合Ⅸ度以上及Ⅴ度以下的经验公式,但是拟合方法与Ⅵ—Ⅷ度的拟合方法是一致的。


图 1 本文仪器烈度算法流程图 Fig. 1 Flow chart of the seismic intensity rapid assessment method
3 计算实例

利用汶川8.0级强震记录,通过本文推荐的仪器烈度算法、袁一凡(1998)提出的算法(以下简称模糊算法)及中国地震烈度表推荐的峰值指标,分别计算仪器烈度值并考查其与实际现场调查烈度的相符性,进而验证该算法的准确性。

金星等(2010)对汶川8.0级地震进行了烈度调查,绘制出烈度分布图(图 2)。其中Ⅵ度以上的地区内强震观测台站有89个,Ⅵ度区内强震观测台站有58个,Ⅶ度区内强震观测台站有20个,Ⅷ度区内强震观测台站有6个,Ⅸ度区内强震观测台站有5个。根据这些强震观测台站强震记录,通过本文推荐的仪器烈度算法进行烈度评定,结果见图 34


图 2 汶川地震现场调查烈度分布 Fig. 2 Intensity isoseismal map of the 2008 Wenchuan earthquake

图 3 3种方法计算结果对比图 Fig. 3 Comparison of calculation results of the three methods

图 4 3种计算方法烈度误差比较 Fig. 4 Error comparisons of the three methods

图 4可以看出,本文算法求得的烈度值与现场调查烈度的相符率为56%,采用模糊算法求得的烈度值与现场调查烈度的相符率为46%,采用中国地震烈度表推荐的峰值指标算法求得的烈度值与现场调查烈度的相符率为45%,本文算法的相符率高于另外2种算法。

宏观烈度如用震害情况评定,评定方法具有多指标性、模糊性、平均性、主观性等特点(胡聿贤,1988),鉴于此,认为仪器烈度值与实际现场调查烈度值相差±1度均属于可接受范围。故从表 2可知,本文推荐的算法具有较高的准确性和可靠性。

表 2 判定结果统计 Table 2 Statistical analysis of the three methods
4 结论

本文通过对大量强震数据统计分析,选取合适的滤波频带,提出了采用三分向合成加速度、速度计算仪器烈度,并利用现有的强震数据给出了Ⅵ—Ⅷ度的拟合烈度公式。同时,利用本文的算法计算得到汶川地震相关台站的仪器烈度,计算结果与现场烈度值具有较高的相符率,证明本文推荐的仪器烈度算法具有一定的可靠性和准确性。但由于用于拟合烈度经验公式的强震数据主要集中在中国西部地区,同时,用于验证本文算法的强震数据也属于中国西部地区,因此,若想得到适用性更广泛的仪器烈度算法,还需要积累更多地区的强震数据,并对拟合公式进行不断的修正、完善。

需要说明的是,本文推荐的仪器烈度算法仅为地震烈度速报提供了一种思路,若在核电厂、高铁地震仪表系统内进行应用,还需要进一步深入研究和大量的试验、分析,如地震的哪些参数与核电厂破坏具有较高的相关性等。

参考文献
陈志高, 黄俊, 2015. 核电站地震监测系统及其开发进展[J]. 大地测量与地球动力学, 35(6): 1065-1068, 1073.
冯德益, 林命周, 吴国有, 等, 1989. 地震烈度的两种模糊评定方法[J]. 地震工程与工程振动, 9(2): 45-56.
何召壮, 2012.从日本福岛核泄漏事故看跨国核污染的国家责任.青岛: 中国海洋大学. http://cdmd.cnki.com.cn/Article/CDMD-10423-1012504255.htm
胡聿贤, 1988. 地震工程学[M]. 北京: 地震出版社.
金星, 张红才, 李军, 等, 2010. 地震仪器烈度标准研究报告[M]. 福州: 福建省地震局.
金星, 张红才, 李军, 等, 2013. 地震仪器烈度标准初步研究[J]. 地球物理学进展, 28(5): 2336-2351.
李亮, 2011. 基于地震动参数的烈度计算方法研究[M]. 哈尔滨: 中国地震局工程力学研究所.
李敏, 2010. 地震动加速度反应谱与地震烈度的关系研究[M]. 哈尔滨: 中国地震局工程力学研究所.
李山有, 金星, 马强, 等, 2004. 地震预警系统与智能应急控制系统研究[J]. 世界地震工程, 20(4): 21-26. DOI:10.3969/j.issn.1007-6069.2004.04.003
刘锡荟, 王孟玫, 汪培庄, 1983. 模糊烈度[J]. 地震工程与工程振动, 3(3): 62-75.
苏经宇, 周锡元, 谭健, 1991. 地震烈度的模糊预测模型[J]. 地震学报, 13(3): 387-394.
袁一凡, 1998. 由地震动三要素确定地震动强度(烈度)的研究[M]. 哈尔滨: 国家地震局工程力学研究所.
郑华, 魏淑虹, 2015. 核电厂地震自动停堆综述[J]. 核科学与工程, 35(4): 664-674. DOI:10.3969/j.issn.0258-0918.2015.04.013
周正华, 周雍年, 赵刚, 2002. 强震近场加速度峰值比和反应谱统计分析[J]. 地震工程与工程振动, 22(3): 15-18. DOI:10.3969/j.issn.1000-1301.2002.03.003
Brackman T. B., Withers M., 2006. Implementing ShakeMap for the New Madrid seismic zone[J]. Seismological Research Letters, 77(4): 445-452. DOI:10.1785/gssrl.77.4.445
Kuwata Y., Takada S., 2002. Instantaneous instrumental seismic intensity and evacuation[J]. Journal of Natural Disaster Science, 24(1): 35-42.
Wald D. J., Quitoriano V., Heaton T. H., et al, 1999. TriNet "ShakeMaps":rapid generation of peak ground motion and intensity maps for earthquakes in Southern California[J]. Earthquake Spectra, 15(3): 537-555. DOI:10.1193/1.1586057