地震应急自动触发计算产出历史地震活动性资料及前兆观测信息的研究

洪旭瑜,林苗禄,关玉梅,李 强

(福建省地震局,福州 350003)

0 引言

随着信息化建设日益成熟,依托计算机软件编程、WebGIS、数据库以及实时融合等技术实现定制需求包括自动产出在很大程度上都能够实现[1-3]。林岩钊等[4]基于地震预警与烈度速报系统实现自动产出更接近现实的地震专题图;张博[5]利用 MATLAB自动产出地震活动性资料;刘坚等[6]采用自编程序,实现地震应急与科学产品后台自动加工处理与发布。因此,地震应急及科技产品能够由软件程序实现的,让人工操作向自动化转变是未来发展趋势。

福建已建成前兆台网技术观测系统,图1显示福建省前兆台网各学科观测台站,到2018年底共有270套左右前兆观测仪器,包括形变、重力、流体、电磁等学科的仪器,仪器观测测项数达500个左右。随着前兆仪器数据量及信息日益增多,这使得日常和应急会商需花费较多人工精力对前兆曲线绘制及质量分析,且图件绘制无统一规格标准,而地震应急会商时要在规定时间内快速给出初步的震情专报及会商意见。因此,除震后趋势判定外,通过统计并计算震中一定范围内相关历史地震活动性资料,按照相对固定格式自动产出初步震情专报,同自动绘制产出的规范化前兆观测数据曲线及异常干扰情况说明相结合,为震后趋势判定提供参考,并可应用于日常会商资料产出和各类图件报告的制作。

图1 福建省前兆台网各学科观测台站Fig.1 Observation stations of all discipline of Fujian precursory network

Java发展到现在已经很成熟,它主要由编程语言 (Java语言)、运行环境 (JVM:Java虚拟机)、框架 (Java API群)组成[7],具有面向对象、可移植性、高性能、多线程、动态性等特点,相比MATLAB,JAVA在可视化界面、兼容性和跨平台性具有较大的优势。MATLAB在算法开发、数据可视化、数值计算等高级技术计算语言和交互式环境拥有强大的功能。因此,自动产出程序设计在面向对象可视化展示及算法程序设计语言分别选择JAVA技术和MATLAB编程语言,充分利用二者之间的优势。

考虑用户使用的便利性,采用JAVA调用MATLAB程序,JAVA作为可视化界面进行设计,通过JVM实现平台无关性（可同时兼容32位及64位 Window7-10系统及各版本 Linux操作系统）,无需担心用户使用的操作系统版本,便于系统的安装使用。操作系统在没有安装 MATLAB相关版本是无法使用内置于JAVA的MATLAB程序,但是让每个用户单独安装较大文件 MATLAB程序显然不太现实,采用将MATLAB程序代码打包成 jar文件内置到Java作为一个包,安装 MATLAB软件提供的精简版MCR（MATLAB Compiler Runtime）软件即可实现JAVA与MATLAB之间的代码共享兼容。通过JAVA可视化界面调用MATLAB程序,便捷、高效产出美观、规范化、标准化的会商资料。

1 测震信息产出及其计算原理

1.1 地震活动性参数

地震发生后,震后趋势判定需要结合地震的震中位置,震中一定范围内历史地震分布情况,历史地震序列类型,序列类型时空分布情况等信息进行综合分析。通过地震速报信息系统（EQIM）正式报与省级地震台网速报内容相结合,选择合适的地震目录,算出震区周边历史中强震活动性情况。

1.2 相关参数计算

震中附近历史地震分布情况：需要产出震中附近50 km和100 km范围内历史地震分布情况：计算震级在3.0～3.9级、4.0～4.9级、5.0～5.9级、6.0～6.9级以及7.0级以上地震个数,要得出50 km和100 km范围内发生的最大地震的时间、地点以及震级。因此需要计算历史地震与当前震中的震中距。震中距有多种算法,不同算法精度不一样,其中精度较高、采用较广的算法主要有Robbins法、Gauss法以及地心纬度法,Robbins法、Gauss法对震中纬度（高纬度或是低纬度）以及方位角的要求较高,而计算800 km以内的震中距,受方位角影响较小的一般采用地心纬度法,获取两个点的经纬度就可得出震中距,公式如下[8]：

式中：lat2为历史地震震中纬度,lat1为当前地震震中纬度;long2为历史地震经度,long1为当前地震经度,D为震中距。

历史地震序列类型和序列类型时空分布情况计算原理同历史地震分布情况,这里不再叙述。

1.3 自动触发

自动触发是地震发生后以EQIM系统产出的正式报震级为准,当达到设定的震级,系统自动触发快速计算预先设定好的相关地震活动性参数和信息,绘制震中附近一定范围内的长期、中期、短期前兆观测数据曲线以及对应的异常变化数据曲线情况说明,然后将计算结果和绘制好的曲线及异常曲线情况说明以PowerPoint、Word等报告形式产出。判断自动触发的要素为地点和震级,二者同时满足条件即可触发。例如当福建省行政区内发生M3.0（这里不区分ML和MS）级及以上地震,自动触发条件满足。借用MATLAB函数库timer函数定时刷新功能定时监测EQIM系统正式报信息,假设LastTime为EQIM最新正式报地震的发震时间,BeforeTimew为本地存储前一次地震发震时间,则监测过程如图2。

t=timer('TimerFcn','m文件','Period',循环间隔,'ExecutionMode','fixedSpacing','TasksTo Execute',inf);%创建自动触发定时器;

start(t);%开始执行

由于前兆仪器观测测项较多,自动绘制长中短期需要一定的时间,因此需要考虑整个程序执行完一次需要耗时多长,根据测试结果需要4.5 min左右,循环间隔要大于运行完成时间,同时兼顾时效性,因此设定为5 min。循环次数需设置为无限次循环inf。

图2 会商资料自动触发监控过程Fig.2 Monitoring process automatically triggered by the consultation information

2 前兆观测信息产出

当前绘制前兆观测曲线并形成 PPT会商报告主要是通过EIS2000（地震前兆信息处理与软件系统）和中国地震前兆台网数据处理系统[9-10]筛选台站和测项绘制曲线,并查找相应时间段的观测曲线数据异常干扰情况说明,该方法需耗费较多的人工精力,在产出图件方面可能存在曲线坐标格式、美观度及标准化等方面不统一。MATLAB在绘图方面的强大功能可以较好的解决该问题。

2.1 数据库基础信息

福建前兆台网基础信息数据库采用的是分布式Oracle数据库,其遵循《地震前兆数据库结构》（DB/T 51-2012）设计规范[11]。MATLAB自动绘制观测数据曲线及获取异常干扰情况说明需要获取Oracle数据库相关观测数据表、数据字典信息及日志表,周克昌等[11]对地震前兆数据库系统及其相关表信息已经做了较为详细的介绍,这里不在赘述。

2.2 数据库访问机制

MATLAB访问Oracle数据库有三种方式：

（1）OLE-DB (Ob ject Lin king and Embedding DataBase)。

（2）ODBC(Open DataBase Connectivity);

（3）JDBC (Java Data Base Connectivity)。

本文采用MATLAB作为后台程序供Java调用,因此采用JDBC访问机制。Java的Ojdbc14.jar文件提供MATLAB桥接Oracle的作用,将该文件拷贝到$HOMER2018(因版本可能会有所差别)javajar oolbox中,然后打开classpath.txt文件添加语句$matlabroot/java/jar/toolbox/ojdbc14.jar。连接语句如下,连接成功状态见图3:conn=database(数据库实例名,用户名,密码,'oracle.jdbc.driver.OracleDriver',['jdbc:oracle:thin:@',IP,':1521:'])。

图3 Matlab成功连接Oracle数据库状态Fig.3 Status of MATLAB successfully connecting to Oracle Database

2.3 曲线自动成图

参数信息：绘制曲线需要台站名称、台站代码、测项代码、采样率、测项分量名称以及测项单位等六个参数信息,形成单独M文件。曲线自动绘制是通过对这些参数信息的循环读取依次获取观测数据自动实现绘制。MATLAB从Oracle获取观测数据需要解决两个问题：分钟值和秒值观测数据字段是CLOB类,无法直接获取;原始观测数据和预处理数据个别时段有出现缺记,数据库中数值缺记是用“NULL”或“NULLALL”字符表示,MATLAB是无法识别除NAN字符外其他表示数值的字符。对应的解决办法是利用MATLAB提供的getSubString子类将oracle.sql.CLOB转换为String类型,用regexprep函数将“NULL”和“NULLALL”全部替换成“NAN”。

曲线成图：自动绘制后的曲线是系统默认模式,需要按我们需求对横坐标和纵坐标刻度值、纵坐标单位以及横坐标标题的标注进行规范,规范后曲线如图4。

（1）横坐标刻度：根据数据类型标注时间格式,原始数据和预处理数据时间精确到日,时间格式为“YYYY-MM-DD”,整点值和日值数据时间精确到月,时间格式统一为“YYYY/MM”,根据figure大小选择自适应模式等间隔标注刻度。

（2）横纵表标题：标题统一为“台站名称-观测测项分量名称（起始日期-结束日期）”,其中日期格式为“YYYYMMDD”。

（3）纵坐标刻度：刻度标注根据figure大小自适应模式自动标注,需要值得注意的是静水位数值表示的是井口距离探头的数据,静水位数值上升说明水位下降,反之水位上升,体应变观测同水位观测原理一致,观测值上升表示向压缩状态转变,观测值越小表示向拉张状态转变。因此,除了体应变观测和静水位测项纵坐标刻度值为反方向,其他测项的纵坐标刻度值离原点越远值越大。

（4）纵坐标单位：MATLAB纵坐标单位默认显示在左侧中间,规范统一为纵坐标正上方,语句为：text(min(xdata),axis(4),unit,'HorizontalAlign ment','center','fontsize',15)。

图4 规范化观测曲线图件Fig.4 Curve diagram of standardized observation

3 结果产出

根据对地震活动性参数的计算结果以及曲线自动绘制成图结果分别产出初步的震情专报和前兆各学科观测曲线及其异常干扰情况说明分析报告。Office应用程序套件具有丰富的功能,调用 Office VBA方法可以快速实现对文档、电子表格和演示文稿进行创作、格式设置和操纵。通过 MATLAB调用 office VBA的对象模型、属性和方向实现对 PPT和 WORD模块的设计,实现报告的自动化产出。

3.1 震情专报产出

初步的震情专报自动产出是通过MATLAB调用本地WORD服务器,以预先设定好的相对固定格式和内容形成的模板产出得到专报内容,假设2019年10月25日福建莆田仙游县发生3.7级地震,模拟自动快速产出初步震情专报,产出结果见图5。

图5 模拟震情专报自动产出结果—以2019年10月25日福建莆田仙游县3.7级地震为例Fig.5 Automatic output results of simulated earthquake situation special report --taking the Xianyou M3.7 earthquake in Putian,Fujian Province on October 25th,2019 as an example

3.2 观测曲线及异常干扰情况说明产出

MATLAB调用PowerPoint服务器,将绘制完成观测曲线及曲线异常干扰质量情况说明自动导入PPT文件,一张幻灯片版式由标题和绘制好的观测曲线组成,标题以“台站名称-仪器名称-短（中或长）期”命名,其中短期为曲线时长小于3个月,数据为原始采样值（在没有特殊要求下默认采用观测仪器原始产出值：秒值、分钟值、整点值、日值,下同。）,中期为时长大于三个月小于6个月,数据一般采用整点值,长期为时长大于6个月以上,数据一般采用日值,备注页说明曲线异常干扰情况,导出结果见图6：

图6 观测曲线及异常干扰情况说明自动导入PPT结果Fig.6 Observation curve and abnormal interference description automatically imported into PPT

4 结语

本文主要利用MATLAB软件在数据处理和计算能力方面的强大功能,能够快速计算相关参数结果及占用较少的内存获取较大字节的地震前兆观测数据,保证部署在服务器端实时监测地震发生后快速判断是否满足条件并自动产出相应可靠的资料,争取应急时间,采用JAVA对程序内容进行可视化,满足日常会商对图件和资料个性化产出的需求,极大的提高工作效率。