基于谱比法的台站场地响应软件研制

汤曜玮，高 亮，邓远立，陈 瑛，龚 萱，曾晓燕

（广东省地震局，广州 510070）

0 引言

一般来说，固定台站均选址在地质条件较好的基岩之上，可以认为台站的场地响应近似取1，但是也有研究结果表明，基岩场地并不是完全地无放大作用。而在台网的加密过程中，为了台网布局的需要，有些台站的场地响应也并不是很好，在大震应急中，也有很多流动台站并不是布设在基岩上，而是在软土层上。对于松散土层，由于其密度和波速相对较低，造成介质的阻抗较低（介质的阻抗等于密度和波速的乘积），而地震波的振幅和阻抗的平方根成反比，因此布设在松散土层上的流动台相比与理想基岩，其地震计振幅会更大，基于此，在地震研究中，考虑对应台站的场地响应，研究场地响应计算方法，开发了与之对应的操作简单、实用化的场地响应计算软件，有助于基层地震工作者用于新建固定台站、流动台场地响应测算或其他工程地震研究方面，具有一定的应用价值。

Matlab 是由MathWorks 公司推出的一款强大的数值计算软件，是数值计算、符号运算和图形处理等多种功能的强有力实现工具[3]。近年来Matlab这一强大的科学计算软件包已经得到业界的广泛认可，并已深入到了各行各业的众多学科、在各大公司、科研机构、大学校园得到了日益普及与广泛应用。其具有以下优点：①代码简捷直观，以矩阵为运算基础，在计算方面避免使用者花费大量时间去编写不必要的小程序；②数值计算方面，Matlab基本涵盖了所有的数值计算方法，从基本的求极限、解微分方程，数字信号处理FFT，到有限差分、线性拟合，再到最新的最优化技术诸如神经网络、遗传算法等等，使复杂的数值计算变得容易；③数据可视化方面，Matlab提供了多种基本图件处理函数，极大的方便绘制出高质量的工程图谱；④Matlab提供了强大的编程接口，支持Matlab与其他编程语言进行数据交换。自地震观测从模拟观测步入数字化观测以来，地震仪器拾取的地面震动信号本质上也是一种随机数字信号，而Matlab编程平台自带的数字信号处理工具箱对于数字信号的各类型常规处理诸如卷积、滤波、加窗、FFT、相关性分析等内置了多种函数[4]，为本软件的研制也提供极大的编程方便，节省了开发时间。

1 谱比法场地响应计算原理及数据处理方法

1.1 谱比法场地响应计算原理

设VS为地表处地震动垂直分量的地震波振幅谱，VB为基底处地震动垂直分量的地震波振幅谱，HS为地表处水平分量的地震波振幅，HB为基底处的水平分量的地震波振幅谱，则经验传递函数[5]为

即井下（基岩处）地震动的水平分量和垂直分量大致相等，垂直分量经过软土层的放大效应远远小于水平分量的放大效应，基本没有被放大，在此基础上，经验传递函数可以简化为[6]谱比率S=SN，于是有

根据上式（6）用地表记录到的水平向地震波振幅谱和垂直向地震波振幅谱比就可以计算得到各台站的场地响应，根据以往文献，除了利用地震事件计算场地响应，同样也可以利用地震计记录到的地脉动来计算场地响应。在本软件中，可以对地震事件和地脉动分别进行计算求取对应台站场地响应。

1.2 数据处理方法

当采用台站记录到的地震事件来计算场地响应时，在对地震记录数据进行谱分析时，由于不同震中距的台站记录到的同一地震事件，其地震记录S 波持续时间并不相同，如果对记录S 波直接进行FFT 转换求取幅度谱，其频率谱间隔疏密必定不同，为了得到相同间隔的频谱，在软件中，我们采用经典平移窗谱[7]的方法。即把所取的S 波窗分割成若干个包含256个数据点的小段，并使相邻数据段有50%的重叠间隔，以获得稳定的波形记录的傅里叶谱，同时为了防止切割波形记录产生频谱泄露，在每一小段波形的起始和末尾加以5%的COS 边瓣。对于采样率为100 Hz 的地震记录来说，每个小段的时间长度是2.56 s，通过快速傅里叶变换得到每个小段的傅里叶谱。然后对每个小段的傅里叶谱进行仪器响应校正，这样通过下式[8]

我们可以得到整个“S 窗”内观测信号的位移谱振幅，式中v2i( f )是经过校正掉仪器响应之后的S波波形记录分割后第i小段的傅里叶快速变换谱，T 为“S 窗”的持续时间，该S 窗包含了n 个时间长度为t 的256 个采样点的小段。最后然后通过内插得到对数频率为0.0， 0.05， 0.1， 0.15， 0.2，…，1.3 的共计27 个频率点的位移谱振幅值。

在以上进行谱分析之前，由于地震台站地震记录中均可能出现直流偏移，这些偏移均可能包含在地震事件记录中或者使地脉动波形记录不是真正的地脉动。在实际计算中均需要进行直流分量剔除。

角砾状构造：早期形成的黄铁矿被错动成棱角分明，大小不等的角砾分布在石英脉中，有的角砾被次生铁锰矿物胶结。

2 软件设计思路及部分实现方法

根据上述谱比法场地响应计算原理以及平移窗谱法数据处理方法，软件总体设计思路可分为以下几个步骤:①导入原始地震事件波形数据或地脉动波形；②导入对应地震台站仪器参数数据；③导入地震事件震相数据，根据P 波与S 波走时差经验估算S波持续长度，对于地脉动记录则设置统一时间长度；④计算某批次地震各台三分量位移振幅谱；⑤合成水平向振幅谱比垂直向得到某批次记录各台场地响应值；⑥对多批次地震记录或地脉记录求得场地响应求均值。

2.1 数据导入

软件导入的地震记录数据为JOPNES 地震速报分析系统截取的导出的sac 数据文件，sac 数据格式广泛用于地震学分析软件中，sac 格式也是地震波形数据的通用格式。sac 数据格式有ASCII 码和二进制两种形式，本软件数据导入采用sac 二进制码格式。软件实现了对多台地震事件记录的批量读取。在批量读取地震事件记录的同时，也批量读入对应地震台站地震计参数信息，其地震计参数信息是JOPNES 地震速报分析系统截取的导出的ascii 文件中的par 文件，其中包含了对应地震记录的地震计灵敏度、传递函数零极点等信息。另外，对于地震事件，为了准确确定地震记录S波流逝事件，软件也直接批量读入JOPNES 地震速报分析系统产出的震相phase文件。

2.2 地震S波窗长确定

精确的S波持续窗长是地震记录位移谱分析的基础，在软件中导入的phase 文件均是对应地震记录不同台站各类波到时的编目结果，精度可靠。在导入phase 文件后，对于S 波持续时间长度的确定，根据刘杰等指出用于计算位移谱的S波时间长度应为为从S 波开始到包含90%S 波能量的时间段[9]，通过回归分析1900 条地震记录得到S 波截取时间长度与Pg、Sg 波走时差为线性关系，并拟合出具体线性表达式，如下式所示：

其中Ttotal是S波截取持续时间，TSg是S波到达时间，TPg是P波到达时间，在软件图形化波形显示区域，根据导入的phase文件编目震相和走时经验公式其S波窗长是可靠的。

2.3 场地响应值求取

对所计算台站根据其每次地震记录三分向利用公式（7）计算地震波位移振幅谱，然后将水平EW 与NS两分向利用公式（9）做合成，垂直向VS（f）即为地震计垂直向位移振幅谱。

再将（9）与（10）求比值，即得该次地震事件不同频点的谱比法场地响应值。对于某个台站，可以计算多次不同地震事件的场地响应值，最后将多批次地震事件场地响应统计求均值即为该台站场地响应。

如果采用地脉动记录来计算场地响应，无需考虑地震事件中S波持续事件问题，我们将平静时段的地脉动数据统一按600 s 一段截取[10]，也是每段数据叠加前一段50%记录，分别FFT 计算每段速度谱并转换位移振幅谱后，每段水平向谱值比垂直向谱得到每段场地响应，最终求得该台站场地响应均值。

2.4 坏谱结果剔除

在计算过程中发现，对于批量处理地震事件时候，部分台站由于漂移、地震记录中断等原因，计算得到的谱值明显异常，为了避免异常谱值影响最终多批次场地响应平均计算结果，在软件中设计了谱值存储前剔除功能。

3 软件使用的简单说明

3.1 软件基本界面

软件主界面如下图1所示，界面分为包括以下几个部分，参数输入及操作流程控制部分，波形显示部分，位移谱值结果及场地响应结果显示三部分，面板左边主要是参数输入及计算流程控制部分，中间三个坐标轴是波形记录三分量显示部分，左边两个坐标轴用于显示位移谱值结果与场地响应结果。参数输入部分又分别按照地震记录或地脉动记录两部不同计算方式进行分类输入。在左边的操作流程控制框里，依次导入各项参数数据后，点击生成位移振幅谱即可在软件后台自动完成位移谱计算全过程，然后通过下拉框可以选择对应台站，在中间及右边观察对应台站三分向波形记录或地脉动记录、位移谱和场地响应结果、根据观察情况选择是否保存该次谱结果，在多批次的结果计算完成之后，可以存档保存，最后点击第三个面板框里的求取场地响应均值完成最终计算并保存结果图件。

3.2 软件操作流程

分别以地震事件和地脉动记录为数据源介绍软件使用流程，软件的基本功能架构图如下图2所示，首先导入各类基础数据，后台计算，结果显示，结果重处理。

图1 程序界面Fig.1 Program interface

图2 基本功能架构图Fig.2 the basic function diagram

当处理地震事件记录时，点击左侧上方面板框内批量导入地震事件、批量导入震相文件、批量导入仪器参数完成后，再点击批量生成场地响应，即可在软件完成该批次地震所有台站的场地响应计算，场地响应包含的如S波窗长确定、FFT谱计算、谱比法计算过程都在软件内部后台自动计算，计算完成后，在右边上部坐标轴显示出该次地震所有台站的水平与垂直位移谱线，同时在按台站显示的下拉菜单中会自动导入该次计算中的所有台站，然后点选某各台站，可以在右中部显示坐标轴中显示该台站此次地震波形、水平及垂直向位移谱结果，如下图3所示。

中间三坐标轴显示德庆台该次地震记录，并用三条红线显示了PG、SG到时与SG波截断时长。右侧下部坐标轴为该次地震记录位移谱。同时软件会将该次地震所有计算台站位移谱结果保存至硬盘某文件夹，文件保存内容如下图4 所示，含有该次地震相关信息，各台的不同频点位移谱值通过软件的显示界面，我们可以观察求得的位移谱，如果发现有某台的谱值异常，我们可以选择剔除坏谱，不予保存至文件里。

当计算多批次地震之后，多批次的地震位移谱结果会依次保存下来，然后点选软件面板左部最下方多批次结果求平均的控制框，读入保存文件里的震幅谱结果，再选择台站既可以得到该台多次记录结果的场地相应的平均值，如下图5 所示，图四上方坐标轴为两次地震记录的震幅谱，其中粗线为其震幅谱均值，下方坐标轴为其谱求平均之后，然后进行频比求得场地相应，其中绿线为最终场地响应。地脉动记录软件操作过程与之类似，在此不在赘述。

图3 德庆台地震波形记录、该次地震位移谱及场地响应Fig.3 Seismic waveform record，displacement spectrum and site response of the earthquake recorded by Deqing station

图4 地震记录位移谱结果文件Fig.4 Displacement spectrum result file of seismic records

图5 震幅位移谱均值及场地响应结果Fig.5 The mean value of the seismic amplitude displacement spectrum and site response results

图6 部分台站场地响应结果Fig.6 Site response results of some stations

4 软件使用实例

选取广东台网2017—2018年间7次发生在省内速报地震，所有地震震级均在ML2.5 级以上，广东台网大部分台站以及部分省外台站均有清晰的记录，按照以上介绍的软件使用流程，依次批量导入这7 次地震的多台地震记录、震相记录及仪器参数，最后计算得到广东大部分测震国家固定台站的场地响应，现选择部分台站结果显示如下图6所示:从结果看出大部分台站场地响应结果在100上下左右浮动，比较符合一般性实际台站场地响应结果预期。说明软件计算场地响应达到预期效果。

5 结论与探讨

本文以Matlab为开发平台，编制了基于谱比法台站场地响应计算软件，该软件界面友好，批量导入地震记录、仪器参数以及震相记录，在导入基础数据后，在软件后台计算各台位移谱，并用文件形式保存结果，最后将某台多次位移谱结果求平均后计算该台场地响应，整个计算过程在软件中集成一体，对于地震工作者来说，该软件界面友好，操作简单，可用于新建固定台站、流动台场地响应测算或其他工程地震研究方面、具有一定的应用价值。