地震前兆信息量计算软件研制及其操作说明1

2013-11-26平建军张永仙单连君平彩鹏

震灾防御技术 2013年4期

平建军张永仙单连君平彩鹏

1）河北省地震局，石家庄 050021

2）中国地震台网中心，北京 100045

引言

目前地震综合分析主要有两大方法：一是从地震前兆异常做起，代表性方法有地震综合计算预报（罗兰格，2002）、综合积分法（李广鑫等，1992）和专家系统（王炜等，1996）等；二是直接从前兆观测数据做起，代表性方法有地震前兆信息量综合数学表达式（平建军等，1999）、综合异常指数（王海涛等，2002）。其中，地震前兆信息量综合数学表达式方法不仅可提取地震异常持续期间的前兆信息，还可提取异常结束后仍延续有效期间的前兆信息，不仅可适用于对各种物理量纲的观测时序数据进行处理，还可使提取到的前兆信息以无量纲方式得到完全统一，从而为开展时间、空间地震综合分析奠定了基础。该方法在提取地震前兆信息时，要求观测数据是平稳时序序列。而实际观测资料中，往往存在着趋势上升、下降、周期等非平稳变化，这就需要选择平滑滤波、傅里叶变换、变化率等数学方法，事先进行消周期、消趋势项的数据预处理工作。我们知道地震前兆信息可能存在于不同的频带中，既可能存在于高频段，也可能存在于低频段，还可能兼而有之。在现行的数据预处理工作中，对前兆观测数据或去低频，或消高频，通常会导致遗漏一些有用的异常信息，使异常信息不能完全表征出来。近几年越来越被重视且得到广泛应用的小波和经验模态（EMD）分解技术（程正兴等，1994；武安绪等，2006），为解决这一问题带来了便利。通过对原始前兆观测资料进行小波或EMD分解，可有效地把各种频率成份从中分离出来，且分离出来的诸时频序列资料数据变化平稳，非常便于应用地震前兆信息量综合数学表达式对其进行无量纲化震兆信息提取。

本研究将地震前兆信息量综合数学表达式和小波、EMD分解技术结合在一起，分析处理前兆观测资料，并据此研制了相应的计算软件，以期为多频域剖析前兆观测资料的变化特征、客观自然地多层次挖掘地震异常、最大程度地展示前兆群体异常映震效能，提供一个实用性、可操作的技术平台。

1 研究思路简介

对工作区域每一前兆观测资料（包括：地壳形变、地电、地磁、水化、水位、应力、应变等）首先判别其是否存在非平稳变化。若存在非平稳变化，则使用小波或经验模态（EMD）分解技术，将其分解为一组多阶时频序列，分析各时频序列映震能力，并挑出映震效能最优的时频序列，应用地震前兆信息量综合数学表达式，提取其无量纲前兆信息量；如果不存在非平稳变化，则可根据地震前兆信息量综合数学表达式，直接提取计算其无量纲前兆信息量。地震前兆信息量综合数学表达式及其式中各参数的含义，请参见有关文献（平建军等，1999），这里不再复述。

对工作区域各种量纲的前兆观测资料，统一按上述方法提取前兆信息量后，即可开展如下时空综合分析。

1.1 地震综合前兆信息量时序分析

根据每一前兆观测资料地震信息量的提取计算结果，按下式逐时刻地截取并计算工作区域地震综合前兆信息量时序值：

式中，S时序为某时刻工作区域地震综合前兆信息量时序值；Si为该时刻某前兆观测资料的前兆信息量值；N为该时刻所有前兆观测资料总数。

按式（1）进行截取计算后，可得到一条工作区域时序综合前兆信息量曲线，据此可研究并提取地震时间综合预测判据指标。

1.2 地震综合前兆信息量空间分析

根据工作区域每一前兆观测资料地震信息量的计算结果及其台站的空间经纬度，按式（2）逐时刻地截取并绘制工作区域空间综合前兆信息量分布图，通过研究强震前区域震兆信息空间演化特征，从而建立其强震空间发震地点判据指标。

式中，S窗为某时刻工作区域某经纬节点扫描范围内截取并计算的地震综合前兆信息量值；Si为该时刻该经纬节点扫描范围内某前兆观测资料的前兆信息量值；N为该时刻该经纬节点扫描范围内前兆观测资料总数。

确定式（2）中经纬节点扫描范围半径的做法是，首先选定欲研究的地震危险性的目标地震震级，根据目标地震震级与地震前兆异常展布范围的统计关系，得到经纬节点扫描计算信息量的范围半径。传统的经验认为：5—5.9级地震的前兆异常展布范围为200km，6—6.9级地震的前兆异常展布范围为300km，7级以上则为500km（罗兰格，2002）。蒋海昆等（2009）由震例统计给出的震级M与最大平均异常震中距D的关系式为：

依据以上所述，可视具体情况，选择相应方法得到经纬节点扫描范围半径。

另外，式（2）中的扫描量S窗也可用经纬节点扫描范围内的异常台项比来表示，即亦可把式（2）改写为

式中，S窗为某时刻工作区域某经纬节点扫描范围内截取并计算的地震综合前兆信息量值；n为该时刻该经纬节点扫描范围内前兆观测资料有前兆信息显示的资料总数；N为该时刻该经纬节点扫描范围内前兆观测资料总数。

2 计算软件设计与操作说明

根据以上研究思路，笔者设计了地震前兆信息量计算软件的数据处理流程（图1），进而研制了地震前兆信息量计算软件。

图1 地震前兆信息量计算数据处理流程Fig. 1 Data management procedures of calculation on earthquake precursory information

该软件是在Windows XP操作环境下，用Visual Basic语言编写的，由13个子窗体、19个计算模块组成，分为工作区域、数据处理、前兆信息量计算、前兆信息量综合分析、显示打印等5个主控菜单，其主控界面如图2所示。下面对5个主控菜单的主要功能、操作说明简述如下：

2.1 工作区域

图2 地震前兆信息量计算软件主控界面Fig. 2 Main control menu of the software

工作区域下属子菜单为选择前兆资料分析工作区域，其相应计算模块为m11。选择前兆资料分析工作区域的功能是：供用户选择输入工作区域的名称及经纬度范围等信息。

启动m11计算模块，首先会出现一个选择或确定工作区域的对话窗，其上列有首都圈地区、川滇地区、青藏高原东北缘和其他地区等4个选项。其中，首都圈地区、川滇地区、青藏高原东北缘是特定工作区域，这些工作区域的经纬度范围数值已直接写入程序，无需用户输入，可直接点击选定。当用户的选择为其他地区时，计算机接连弹出2个对话窗，分别让用户输入工作区域的名称以及经纬度范围。

一旦选定新的工作区域，计算机还需要用户回答，是否导出、删除前工作区域信息量计算的数据处理、计算结果等有关信息，是否导入本工作区域以往信息量计算的数据处理、计算结果等有关信息。

完成以上操作后，计算机即在本软件应用程序目录下创建Gzqyfile子文件夹，在该子文件夹中，形成存储工作区域名称编号、工作区域经纬度范围数值的 gzqy.txt文件；如果用户选择的是其他地区，则在Gzqyfile子文件夹中，还会形成存储工作区域具体名称的diming.txt文件，供程序有关计算模块调用。

2.2 数据处理

数据处理包括数据准备、小波分解和EMD分解3个子菜单，其相应计算模块分别为：m21、m22、m23。

（1）数据准备

数据准备的功能是：调取工作区域单项前兆观测资料，并选定采用其哪段数据作为正式应用数据文件，进行地震前兆信息量计算。

启动 m21计算模块，首先出现一个选择输入数据文件对话窗，此时可通过选择驱动器、文件夹，然后点击有关数据文件名，即方便地把欲选取的原始前兆观测资料数据文件名输入计算机。这里需指出的是，原始前兆观测资料数据格式均为：日期、数据。日期用年、月、日8位数字表示，如某观测资料1989年5月的月均值数据为34.5，则其数据格式为：19890531 34.5。

接着计算机出现更改数据文件名窗口，在其上显示两个文本框，一个是“数据文件名为”文本框，另一个是“数据文件更名为”文本框，这两个文本框的缺省值均为选取输入的原始前兆观测资料数据文件名。本软件要求数据文件名由 8位以内的英文字符或数字加后缀 txt组成，当不符合这一要求时，需对其进行更名。

最后计算机弹出一个对话窗，对话窗上的开始时间、结束时间，是让用户选择将原始前兆观测资料的哪段数据转为正式应用数据文件，用于进行前兆信息量计算，缺省值为其全部时段；另外对话窗上还需用户回答，在转为正式应用数据文件过程中，是否需要转变原数据类型，如从日均值转为5日均值、旬均值、月均值等。

完成以上操作后，计算机即把原始前兆观测资料数据文件，转化为计算地震前兆信息量的正式应用数据文件，并在本软件应用程序目录下创建Userfile、Tempfile子文件夹，在Userfile子文件夹中形成正式应用数据文件*.txt；在Tempfile子文件夹中形成存储包括原始前兆观测资料数据文件全名（含路径）、正式应用数据文件全名等内容的文件Userfile.nam，供有关计算模块调用。

（2）小波分解

小波分解的功能是：将非平稳变化的前兆观测资料，分解为一组多阶时频序列，分析各时频序列映震能力，并挑出映震效能最优的时频序列，进行地震前兆信息量计算。

启动m22计算模块，首先出现一个选择输入数据文件对话窗，其缺省值即为Userfile子文件夹中的正式应用数据文件*.txt，一般直接确认即可。

接下来计算机让用户回答，欲对已选定的数据文件的哪段数据进行小波分解，缺省值为其全部时段，一般直接确认即可。

然后计算机即把小波分解的各阶时频序列数据曲线展示在屏幕上，并让用户挑选确定映震效能最优的时频序列，用以计算地震前兆信息量。

完成以上操作后，计算机在本软件应用程序目录下创建Xbbhfile子文件夹，在Xbbhfile子文件夹中保存小波分解最优时频序列数据文件*.txt，供提取地震前兆信息量的有关计算模块调用。

（3）EMD分解

EMD分解的功能及其计算模块的操作运行过程与小波分解几乎完全相同，这里不再赘述。需要指出的是，m23计算模块运行完后，计算机在本软件应用程序目录下创建 Htbhfile子文件夹，在Htbhfile子文件夹中保存EMD分解最优时频序列数据文件*.txt，供提取地震前兆信息量的有关计算模块调用。

2.3 前兆信息量计算

前兆信息量计算包括求取均值和均方差、提取计算前兆信息量2个子菜单，其相应计算模块分别为：m31、m32。

（1）求取均值和均方差

求取均值和均方差的功能是：对选定的数据文件的时序数据计算其均值和均方差，为下一步前兆信息量的计算做准备。

启动m31计算模块，会出现一个选择输入数据文件对话窗，缺省值为小波或EMD分解最优时频序列数据文件*.txt、或为正式应用数据文件*.txt，一般回车确认即可。

接下来计算机让用户输入开始时间、结束时间、误差系数等信息，其目的是让用户在已选定的数据文件中，选取数据变化较平稳的正常时段来计算均值和均方差。正常时段的开始时间、结束时间的缺省值，为已选定的数据文件的起止时间；误差系数是指去掉正常时段数据中超过几倍均方差的不规则数据后，再求均值和均方差，误差系数参数的缺省值为9.9。

该计算模块运行完后，将弹出一个计算结果信息窗，同时将在本软件应用程序目录下创建Averfile子文件夹，在Averfile子文件夹中形成存储所计算的均值、均方差等参数的数据文件*.par，供提取地震前兆信息量等有关计算模块调用或计算参数查询。

（2）提取计算前兆信息量

提取计算前兆信息量的功能是：对选定的数据文件的时序数据转换计算其地震前兆信息量，为下一步前兆信息量时空综合分析做准备。

启动m32计算模块，会出现一个选择输入数据文件对话窗，缺省值为求取均值和均方差时选定的数据文件，一般回车确认即可。

接下来计算机让用户输入开始时间、结束时间、误差系数、信息量延续系数、数据类型、台站经纬度等信息。开始与结束时间是指用户欲对已选定的数据文件中，哪段数据进行前兆信息量计算，缺省值为其全部时段；误差系数是指用户欲取几倍均方差作为异常判别标准来进行前兆信息量计算，其缺省值为2；信息量延长系数的缺省值为3；数据类型是指已选定的数据文件的数据是日均值还是 5日均值、旬均值、月均值等，当数据为日均值时输入 1，5日均值时输入5，依次类推；台站经纬度的参数需以度、分的格式输入。

该计算模块运行完后，将在本软件应用程序目录下创建Compfile和Lolafile子文件夹，在Compfile子文件夹中形成所计算的地震前兆信息量数据文件*.res，形成存储有关计算参数的数据文件*.par；在Lolafile子文件夹中形成存储所输入的台站经纬度的数据文件*.spa，供地震前兆信息量综合时空分析等有关计算模块调用。

2.4 前兆信息量综合分析

前兆信息量综合分析包括截取计算工作区域时序综合前兆信息量、制作工作区域时序综合前兆信息量曲线、显示工作区域时序综合前兆信息量数据、截取计算工作区域空间综合前兆信息量、制作工作区域空间综合前兆信息量分布图、显示工作区域空间综合前兆信息量数据6个子菜单，其相应计算模块分别为：m41、m42、m43、m44、m45、m46。

（1）截取计算工作区域时序综合前兆信息量

截取计算工作区域时序综合前兆信息量的功能是：根据式（1），逐时段截取并计算工作区域前兆信息量综合值，为绘制工作区域时序综合前兆信息量曲线做准备。

启动m41计算模块，出现一个挑选数据文件对话窗，计算机自动将本软件应用程序目录下Compfile子文件夹保存的所有地震前兆信息量计算结果数据文件*.res全部开列出来，此时用户可对开列的数据文件*.res，进行删除与添加等筛选操作。

接下来计算机让用户回答，截取并计算工作区域时序综合前兆信息量的开始时间、结束时间以及每隔多长时间（即时间间隔）截取计算一次前兆信息量综合值，回答完提问，计算机即进行截取运算。

该计算模块运行完后，将在本软件应用程序目录下创建Timeresu子文件夹，在Timeresu子文件夹中形成工作区域时序综合前兆信息量结果数据文件timeresu.txt，形成存储有关计算参数的数据文件timeresu.par，供绘制工作区域时序综合前兆信息量曲线调用，或供其他计算模块调用显示、查询相关信息。

（2）制作工作区域时序综合前兆信息量曲线

制作工作区域时序综合前兆信息量曲线的功能是：绘制显示工作区域时序综合前兆信息量曲线，据此可研究地震群体异常前兆信息时序变化特征，提取地震综合时间预测判据指标。

启动m42计算模块，出现一个选择显示时序曲线的对话窗，其缺省值为时序综合信息曲线，直接确认即可。

接着计算机让用户选择输入时序综合前兆信息量数据文件，其缺省值为本软件应用程序目录下Timeresu子文件夹中的timeresu.txt数据文件，直接确认即可。

然后计算机自动搜集有关信息，自动给出曲线缺省图名，用户可根据情况进行修改。

随后计算机让用户输入在图上标注地震的目录文件、标注地震的最低震级、开始时间、结束时间等信息。这里标注地震的目录文件为 wkf格式，开始时间与结束时间是指欲显示timeresu.txt数据文件的哪段数据曲线，缺省值为其全部时段。

最后计算机弹出一个选择并确定绘图参数对话窗，让用户回答纵轴欲标注的最大数值、最小数值、刻度等份，缺省值分别是timeresu.txt数据文件的最大值、最小值及5等份，用户可根据情况调整这些参数。

该计算模块运行完后，将在计算机屏幕上显示工作区域时序综合前兆信息量曲线，并在本软件应用程序目录下创建Curvfile子文件夹，此时点击图上的存图按钮，即可把曲线图形保存到Curvfile子文件夹中，图名是timeresu.bmp。

（3）显示工作区域时序综合前兆信息量数据

显示工作区域时序综合前兆信息量数据的功能是：在计算机屏幕上，显示本软件应用程序目录下Timeresu子文件夹中的timeresu.txt文件数据值，以便用户查询有关计算结果信息。

启动m43计算模块，出现一个输入参数对话窗，对话窗上的开始时间、结束时间，是指欲显示timeresu.txt数据文件的哪段数据值，缺省值为其全部时段。回答完毕，计算机即将数据显示在屏幕上，当数据较多时，用户可操作鼠标进行滚动显示。

（4）截取计算工作区域空间综合前兆信息量

截取计算工作区域空间综合前兆信息量的功能是：截取工作区域内某时间截面上的地震前兆信息量空间分布参数，包括每一前兆观测资料的台站经纬度、信息量值。

启动m44计算模块后，其操作过程与m41计算模块完全相同，这里不再复述。

该计算模块运行完后，将在本软件应用程序目录下创建 Resufile子文件夹，在 Resufile子文件夹中，对每一时间截面所截取得到的地震前兆信息量空间分布参数，以数据文件形式予以保存，每一时间截面的数据文件名由该时间截面的日期（8位数表示）加后缀dat组成，供绘制工作区域空间综合前兆信息量分布图调用。

（5）制作工作区域空间综合前兆信息量分布图

制作工作区域空间综合前兆信息量分布图的功能是：根据已截取的工作区域内某时间截面地震前兆信息量空间分布参数数据文件*.dat，按式（2）计算工作区域每一扫描经纬节点的空间综合前兆信息量值，进而绘制该时间截面工作区域空间综合前兆信息量等值线图。

启动m45计算模块，出现一个输入参数对话窗，对话窗上的预期地震震级，是指欲绘制工作区域未来几级地震的空间综合前兆信息量分布图，震级的缺省值为5.0级。

接下来计算机弹出扫描方式对话窗，其上列有选择经纬节点扫描半径的3种方法，还列有选择经纬节点扫描计算信息量值的2种方式，用户自己可酌情选择。

随后计算机让用户输入扫描计算的区域范围及经纬节点的扫描步长，这些参数均以度为单位，其中区域范围经纬度缺省值，为m11计算模块已给定的工作区域经纬度范围值。

然后计算机调取已截取的工作区域内某时间截面地震前兆信息量空间分布参数数据文件*.dat，对工作区域每一经纬节点的空间综合前兆信息量进行计算。计算完毕，给出一个勾画等值线分布图的绘图对话窗，让用户回答勾画等值线的最大值、最小值、等值线间隔值，这些参数的缺省值已从计算结果中索取给出，用户可根据情况调整这些参数。

该计算模块运行完后，将在计算机屏幕上显示该时间截面工作区域空间综合前兆信息量分布图，并在本软件应用程序目录下创建 Surffile子文件夹，此时点击图上的存图按钮，即可把该图形保存到 Surffile子文件夹中，图形文件名由该时间截面日期（8位数表示，如19971231）加后缀bmp组成。另外，在Surffile子文件夹中，还形成该时间截面工作区域空间综合前兆信息量等值线数据文件spacjg.dat，据此可用其他专业绘图软件，绘制出更加精美的等值线图。

（6）显示工作区域空间综合前兆信息量数据

显示工作区域空间综合前兆信息量数据的功能是：在计算机屏幕上，显示本软件应用程序目录下Surffile子文件夹中的spacjg.dat文件数据值，以便用户查询有关计算结果信息。

启动m46计算模块，即将数据显示在屏幕上，每行显示的内容为经纬度、信息量值，当数据较多时，用户可操作鼠标进行滚动显示。

2.5 显示打印

显示打印包括显示原始数据、显示过程处理数据、显示前兆信息量计算结果数据、打印原始数据、打印过程处理数据、打印前兆信息量计算结果数据、显示原始—过程处理—信息量数据曲线图形7个子菜单，其相应计算模块分别为：m51、m52、m53、m54、m55、m56、m57。

（1）显示原始数据

启动m51计算模块，首先出现一个选择输入数据文件对话窗，其缺省值即为Userfile子文件夹中，由原始前兆观测资料所形成的正式应用数据文件*.txt，一般直接确认即可。

接下来出现一个输入参数对话窗，对话窗上的开始时间、结束时间，是指欲对已选定的数据文件的哪段数据值进行显示，缺省值为其全部时段。回答完毕，计算机即将数据显示在屏幕上，当数据较多时，用户可操作鼠标进行滚动显示。

（2）显示过程处理数据

启动m52计算模块，用户可选择显示本软件应用程序目录下Xbbhfile或Htbhfile子文件夹中的最优时频序列数据文件*.txt的数据值，操作过程与m51计算模块相同，这里不再复述。

（3）显示前兆信息量计算结果数据

启动m53计算模块，显示本软件应用程序目录下Compfile子文件夹中的地震前兆信息量数据文件*.res的数据值，操作过程与m51计算模块相同，这里不再复述。

（4）打印原始数据、打印过程处理数据、打印前兆信息量计算结果数据的功能及相应计算模块（m54、m55、m56）的操作过程，分别与显示原始数据、显示过程处理数据、显示前兆信息量计算结果数据基本相同，不同之处在于，前者把有关数据信息打印在纸张上，而后者则是显示在计算机屏幕上，不再复述。

（5）显示原始—过程处理—信息量数据曲线图形

显示原始—过程处理—信息量数据曲线图形的功能是：将由某原始前兆观测资料数据文件所形成并保存在本软件应用程序目录下，Userfile子文件夹中的正式应用数据文件*.txt、Xbbhfile或Htbhfile子文件夹中的最优时频序列数据文件*.txt、Compfile子文件夹中的地震前兆信息量数据文件*.res，以时序数据曲线图形的方式进行显示，供用户分析对该原始前兆观测资料提取的前兆信息量是否合理，如果不合理则重新返回有关计算环节进行调整计算。

启动m57计算模块，其操作过程与m42计算模块基本相同，不再复述。

3 前兆信息量算例

作为例子，应用地震前兆信息量计算软件，分别对河北昌黎水位和兴济地电阻率NE向日均值观测资料，进行了地震前兆信息量提取与计算。

图3a是1988年7月1日至2010年8月8日昌黎水位原始观测资料日均值数据曲线，对其进行小波分解计算，得到的最优时频序列（第5阶时频序列）数据曲线如图3b所示，计算图3b曲线数据的均值和均方差，得到的=0.0041、σ=0.0073，取h=5，按照地震前兆信息量综合数学表达式，提取计算其前兆信息量，得到的结果数据曲线如图3c所示。由图3c可知，在1995年10月6日河北唐山5.0级地震前，昌黎水位从1995年9月27日开始，出现了显著的短临震兆信息。

图3 昌黎水位地震前兆信息量的计算Fig. 3 Calculation of well-water level from Changli station as earthquake precursory information

图4a是1983年1月1日至2011年9月25日兴济地电阻率NE向原始观测资料日均值数据曲线，对其进行经验模态（EMD）分解计算，得到的最优时频序列（第7阶时频序列）数据曲线如图4b所示，计算图4b曲线数据的均值和均方差，得到的=0.0043、σ=0.0083，取 h=3.5，按照地震前兆信息量综合数学表达式，提取计算的前兆信息量结果数据曲线如图4c所示。由图4c可知，在1983年11月7日山东菏泽5.9级和1989年10月19日山西大同6.1级地震前，兴济地电阻率NE向均有较好的短期震兆信息显示，其中1989年10月19日大同6.1级的震兆信息，得到了相关研究结果的佐证（汪志亮等，1990）。