CN115047516B

CN115047516B - 基于Python的长周期谱烈度计算方法

Info

Publication number: CN115047516B
Application number: CN202210497758.5A
Authority: CN
Inventors: 姚新强; 李通; 曹井泉; 王东明; 林逸; 刘红艳; 纪静; 安立强
Original assignee: Tianjin Earthquake Agency
Current assignee: Tianjin Earthquake Agency
Priority date: 2022-05-09
Filing date: 2022-05-09
Publication date: 2024-02-02
Anticipated expiration: 2042-05-09
Also published as: CN115047516A

Abstract

本发明公开了一种基于Python的长周期谱烈度计算方法。本发明包括下述三个步骤：(1)地震加速度数据预处理：对震后各台站监测获得的地震动加速度数据进行处理；(2)调整谱烈度函数公式上下限时间范围为0.1s～10.1s，计算该范围内的相对速度反应谱的积分值；(3)根据谱烈度不同范围分档表，计算出长周期谱烈度值。本发明提出的谱烈度长周期计算方法基于不同台站的强震记录可以迅速计算台站附近的烈度值，不同台站得出的烈度值以点数据的方式进行存储，再利用普通克里金插值法对大量台站烈度点数据进行插值处理，即可绘制谱烈度分布图，可清晰明了的查看极震区的烈度分布情况，为接下来人员伤亡和经济损失评估做好数据准备。

Description

基于Python的长周期谱烈度计算方法

技术领域

本发明涉及地震灾害评估领域，更具体的说，是涉及一种基于 Python的长周期谱烈度计算方法。

背景技术

Housner谱烈度方法：谱烈度值是指某一阻尼比的单质点体系相对速度反应谱在某一周期区间内时域的积分值，是一个与地表地震动反应谱直接相关的物理量，且强调了地震动与结构自振频率相近部分成份的作用，因而谱烈度SI_0.2与地震烈度相关性较好，能够客观地反映出地震动反应谱和结构破坏程度之间的关系。谱烈度的概念由Housner提出，其含义是在0.1s～2.5s范围内的相对速度反应谱的积分值^[1]：

式中，S_v为相对速度反应谱，T为周期，阻尼ξ一般取0、5％或者20％。

谱烈度算法存在的主要问题及不足如下：

研究发现，谱烈度算法相较于传统仪器烈度算法其结果与宏观烈度值误差率最小，具有很高的可靠性，其中Housner提出的相对速度反应谱作为地震动参数在震后进行灾害评估时也显著优于其他地震动参数。

虽然谱烈度作为可靠的地震动强度表征参数获得了广泛的应用，但是该算法同样存在一些缺点，例如给出的烈度值判定范围较小，只能从Ⅴ度烈度值开始判定且只给出Ⅴ、Ⅵ、Ⅶ、Ⅷ、Ⅸ、Ⅹ六种烈度值的判定，而无法像传统仪器烈度方法那样计算出3.2度、4.7度、5.8度等范围更大、更精细化的烈度值；同时，Housner谱烈度公式积分上下限范围为0.1s～2.5s，目的是为了匹配美国当时实际建筑结构的振动周期范围，该范围周期大约涵盖了1层到20层高的建筑结构，而中国建筑结构形式与美国不相同，不同建筑结构周期范围不同，因而无法直接采用美国建筑结构的振动周期范围。

发明内容

基于上述谱烈度算法存在的问题和不足，本发明综合仪器烈度和谱烈度计算方法的优点，提出了一种基于Python的长周期谱烈度计算方法。

本发明一种基于Python的长周期谱烈度计算方法，包括下述三个步骤：

(1)地震加速度数据预处理：对震后各台站监测获得的地震动加速度数据进行处理；

首先选择为有完整正交三分向：东西分向、南北分向、垂直分向的数据；随后对数据进行基线校正和数字滤波处理；基线校正方法为一般采用记录时间过程减去地震事件前10s内记录地算术平均值方法；数字滤波方法则对加速度每个分向采用数字滤波器进行0.1Hz～10Hz带通滤波处理；

(2)反应谱积分值计算：

考虑到国内多种不同建筑结构类型对应多种不同的振动周期，首先调整谱烈度函数公式上下限时间范围为0.1s～10.1s，计算该范围内的相对速度反应谱的积分值，计算公式如下所示：

式中，S_v(ξ,T)是阻尼比为ξ时的单自由度系统的相对速度反应谱，T为周期，ξ取20％，周期范围取0.1s～10.1s；建筑结构破坏主要是受水平向地震作用的影响，计算谱烈度值SI_0.2时仅取水平两方向上对应谱烈度的最大值为该台站最终的谱烈度值；

(3)根据谱烈度不同范围分档表，计算出长周期谱烈度值：

给出谱烈度不同范围的分档表，根据lgSI所属的范围进行计算，表中I_I为中国仪器烈度表给出的不同宏观烈度下对应的仪器烈度值范围；

表

<Ⅵ

Ⅵ

Ⅶ

Ⅷ

Ⅸ

Ⅹ

Ⅺ

Ⅻ

SI

<15.524

15.524～28.379

28.379～66.527

66.527～210.863

210.863～609.537

609.537～1020.939

1020.939～1396.368

>1396.368

lgSI

<1.191

1.191～1.453

1.453～1.823

1.823～2.324

2.324～2.785

2.785～3.009

3.009～3.145

>3.145

I_I

<5.5

5.5～6.5

6.5～7.5

7.5～8.5

8.5～9.5

9.5～10.5

10.5～11.5

>11.5

同时，本发明根据不同地震烈度所对应的仪器烈度I_I的范围，给出lgSI与改进谱烈度I_S之间的一次函数关系分档公式，其中小于Ⅵ度的公式通过Ⅵ～Ⅺ度六组线性公式计算得出的大量lgSI点数据以一次函数的形式拟合而来；由于目前我国规定最大的宏观烈度为Ⅻ度，因此本发明将lgSI>3.145(即反应谱积分值SI>1396.3)的外延部分对应的改进谱烈度值I_S划定为12度；同时，当SI≤1时，由于反应谱积分值特别小，因此对应的震级及烈度也较小，因此本发明给出的谱烈度计算公式基本要求为SI>1。

式中，I_S为长周期谱烈度值。

本发明方法具有三点优势：

1.长周期。考虑到了国内多种不同建筑结构类型对应多种不同的振动周期范围，同时也尽可能充分利用所获得的地震动加速度数据，新的计算方法将相对速度反应谱的积分周期范围确定为0.1s～10.1s。

2.连续性。使用集集地震、汶川地震和芦山地震三次7.0级以上的地震数据进行拟合分析，给出了线性分档式的谱烈度计算公式，计算的公式具有连续性，而非像之前一样只给出Ⅴ、Ⅵ、Ⅶ、Ⅷ、Ⅸ、Ⅹ度单点式的烈度值。

3.精细化。计算出的谱烈度值结果精确到小数点后一位，结果更加精确，而非像之前一样只给出Ⅴ、Ⅵ、Ⅶ、Ⅷ、Ⅸ、Ⅹ度等颗粒度较粗的烈度值。

本发明提出的谱烈度长周期计算方法基于不同台站的强震记录可以迅速计算台站附近的烈度值，不同台站得出的烈度值以点数据的方式进行存储，再利用普通克里金插值法对大量台站烈度点数据进行插值处理，即可绘制谱烈度分布图，可清晰明了的查看极震区的烈度分布情况，为接下来人员伤亡和经济损失评估做好数据准备。

附图说明

图1是本发明的流程框图

(图中PGA`为经过预处理后的峰值加速度)；

图2是CSMNC提供的地震动数据集(以芦山地震为例)；

图3是每个加速度数据的内容及格式(以芦山地震为例)；

图4是对加速度数据进行基线校正处理

图5是对加速度数据进行数字滤波处理

图6是读取表头信息并保存预处理后的加速度数据至Excle文件中；

图7是预处理后的加速度数据；

图8是预处理后三分向数据进行合并；

图9是计算相对速度反应谱的函数response_spectra()；

图10是精细积分法函数；

图11是根据lgSI的值分档计算得到最终谱烈度值I_S；

图12是根据设置的积分区间以及阻尼比计算SI、lgSI、I_S；

图13是每个台站计算结果SI、lgSI、I_S统计；

图14是鲁甸地震台站分布图；

图15是两种烈度对比图；

图16是烈度差分布图。

具体实施方式

下面结合附图对本发明做进一步说明。

如图1所示，本发明一种基于Python的长周期谱烈度计算方法，包括下述三个步骤：(1)地震加速度数据预处理：对震后各台站监测获得的地震动加速度数据进行处理；(2)反应谱积分值计算；(3)根据谱烈度不同范围分档表，计算出长周期谱烈度值。

实施例1：

1.选择地震动三分向加速度数据集，使用国家强震动台网中心 CSMNC提供的数据格式，如图2和图3所示；

2.使用Python对数据集进行预处理，首先采用记录时间过程减去地震事件前10s内记录地算术平均值方法对其进行基线校正处理，如图4所示；然后对加速度每个分向采用数字滤波器进行0.1Hz～10Hz 带通滤波处理，如图5所示。

创建Excle文件夹，读取每个加速度记录表头信息中的台站代码、经度、纬度、加速度方向等信息，如图6所示，并同时将预处理后的加速度数据保存至Excle文件中，如图7所示。

3.将每所台站对应的三个方向的地震动记录合成到一个excle表格中，如图8所示，同时将表头信息也拷贝过来。在合成过程中需要根据方向是否发成冲突以及台站代码是否一致来综合判断是否台站是否符合有三个方向的记录的要求，如果不符合，则去除该台站记录。

4.三个方向的数据合并完成后，选择EW方向和NS方向的数据进行谱烈度的计算，设置积分区间为0.1s～10.1s，同时设置阻尼比为20％，使用精细积分方法求得相对速度反应谱，如图9和图10所示，对根据设置得时间和阻尼参数求得对应得积分值SI_i。

5.根据每所台站的积分值SI_i计算对应的lg SI_i，并根据烈度分档表计算最终对应的谱烈度值I_S，如图11-13所示。

具体例验证：

一、鲁甸地震基本信息

2014年8月3日云南鲁甸县(27.110°N，103.330°E)发生6.5 级地震，震源深度10km，震中烈度Ⅸ度。本次地震震中位置位于鲁甸县龙头山镇，等震线呈椭圆形，长轴走向为北北西向，根据中国地震局发布的烈度图分析可知，Ⅵ度区及其以上总面积约为10350km²，其中Ⅵ度区总面积约为8390km²，Ⅶ度区总面积约为1580km²，Ⅷ度区总面积约为290km²，Ⅸ度区总面积约为 90km²。鲁甸地震触发监测台站较多，共触发79所台站，收集三分向地震动记录237条，台站之间间距较小且采样数据较为理想，综上所述，本发明选用鲁甸地震作为谱烈度分布图的验证案例，各台站分布位置如图14所示。

二、谱烈度分布图对比分析

本发明对地震数据进行编程计算，得出了79所台站阻尼比为 20％下的谱烈度SI及lgSI，同时使用普通克里金插值法处理，绘制出鲁甸地震的谱烈度分布图，并与中国地震局发布的云南鲁甸6.5级地震烈度图进行详细对比，如图15-16所示。

表3.2 鲁甸地震79所台站谱烈度值

本发明将谱烈度分布图与宏观烈度图叠加对比分析，鲁甸地震最高谱烈度值为7.3度，谱烈度等值线长轴呈NW-SE走向，与中国地震局发布的烈度图走势对比一致；震中位置位于5.7度～6.0度等值线附近，普通克里金插值计算震中位置最高为7.3度，与实地震中评估烈度Ⅷ度相差<1度；79所监测台站中有44所台站分布于谱烈度等值线5.4度以上范围内，10所分布于Ⅵ度圈内，其中2所位于5.7度～6.0度谱烈度等值线内，偏差<1度，5所位于6.0度～7.3度谱烈度等值线内，偏差均<1度；有1所台站分布于Ⅸ度圈内，其所在的谱烈度等值线为6.0～7.3度，偏差>1度。综上所述，谱烈度值与宏观烈度吻合率(偏差±1度以内)约为90％，其结果基本能够反应鲁甸地震灾情分布的整体趋势大致范围，具有一定的科学性和准确性。

Claims

1.一种基于Python的长周期谱烈度计算方法，其特征是，按照下述三个步骤进行：

(1)地震加速度数据预处理：对震后各台站监测获得的地震加速度数据进行处理；

首先选择有完整正交三分向：东西分向、南北分向、垂直分向的数据；随后对数据进行基线校正和数字滤波处理；基线校正方法为采用记录时间过程减去地震事件前10s内记录的算术平均值方法；数字滤波方法为对加速度每个分向采用数字滤波器进行0.1Hz～10Hz带通滤波处理；

(2)反应谱积分值计算：

式中，S_v(ξ,T)是阻尼比为ξ时的单自由度系统的相对速度反应谱，T为周期，ξ取0.2，周期范围取0.1s～10.1s；建筑结构破坏主要是受水平向地震作用的影响，计算谱烈度值SI_0.2时仅取水平两方向上对应谱烈度的最大值为该台站最终的谱烈度值；

(3)根据谱烈度不同范围分档表，计算出长周期谱烈度值：

根据不同地震烈度所对应的仪器烈度I_I的范围，给出lgSI与改进谱烈度I_S之间的一次函数关系分档公式，其中小于Ⅵ度的公式通过Ⅵ～Ⅺ度六组线性公式计算得出的大量lgSI点数据以一次函数的形式拟合而来；我国规定最大的宏观烈度为Ⅻ度，因此将lgSI>3.145的外延部分对应的改进谱烈度值I_S划定为12度；同时，当SI≤1时，由于反应谱积分值特别小，因此对应的震级及烈度也较小，因此给出的谱烈度计算公式基本要求为SI>1；

式中，I_S为长周期谱烈度值。