Matlab在天津市强震动台网烈度速报中的应用①

刘双庆，聂永安，高武平，李雅静

（天津市地震局，天津 300201）

0 引言

烈度速报是一项非常重要的地震灾害快速评估工作，快速产出准确的烈度评估图可以有效地指导地震救援，从而极大地减轻地震造成的二次灾害。随着我国强震动台网的建设与完善，烈度速报工作已获得了长足的进展，并在近几年向地震预警工作方面进行开拓［1］。“九五”期间，中国地震局建立了我国第一个由72个实时数据传输强震台和80个电话拨号数据传输强震台构成的地震动强度（烈度）速报台网。在“十五”期间，中国地震局除了继续对北京市、天津市的地震烈度速报台网进行加密建设外，还在昆明、兰州、乌鲁木齐等城市新建地震烈度速报台网，部分省市地震局也正在计划建设具有地震烈度速报功能的强震台网［2］。烈度速报的形式也从过去的宏观调查逐步向仪器烈度侧重［3－4］，发挥仪器烈度在震害快速评估中的作用。烈度等震线也加入了场地因素并对不同烈度段的PGA、PGV、PGD进行了融合，这些技术方法在美国的Shakemap和日本的UrEDAS都有体现［5－7］。泽仁志玛还就等震线的显示问题进行了讨论［8］，刘双庆［9］则利用蒙特卡罗法分析了强震台网布局对等震线的绘制精度的影响。

随着强震动台网的高密度建设，不同厂家不同型号的强震仪器得到了使用，从而产生了多种数据回收方式及数据保存格式。如何快速地处理这些不同型号不同存储格式的大批量强震数据成为了快速烈度速报的阻碍，就目前天津市强震动台网110个强震台而言，就涵盖了四种类型的仪器。其部分回传的数据还存在着坏死而无法用厂家提供的软件打开的问题，且直接利用厂家提供的软件无法批量处理强震数据并快速归算烈度值。因此本文拟利用已成为国际控制界的标准计算软件Matlab来解决上述的主要难题，并对其中的核心处理技巧进行说明。

1 Matlab及天津市强震动台网简介

1．1 Matlab简介

20世纪70年代，美国新墨西哥大学Cleve Moler为了减轻学生编程的负担，用FORTRAN编写了最早的Matlab。1984年由Little、Moler、Steve Bangert合作成立的MathWorks公司正式把Matlab推向市场。到20世纪90年代，Matlab已成为国际控制界的标准计算软件，并在2011年推出了Matlab2011b（7．13版本）。从 Matlab6．5开始，其主要代码已改由Java语言编写，并与Fortran，C，C＋＋，Java，Vb等语言兼容。由于 MathWorks公司允许世界各地的专家将自己开发的Matlab子程序上传到该公司网站，并由MathWorks公司的工程师对程序代码进行复核，通过全面测试后在新版的Matlab中追加上去，因此目前的Matlab涵盖了多个模块的应用程序，主要包括：Web服务、信息通信、财政金融、控制系统设计、曲线拟合、数据库使用、滤波器设计、模糊逻辑分析、仿生算法、图像处理、地图工具、神经网络、Mu－分析、优化分析、Robust控制、信号处理、样条工具箱、统计工具箱、符号数学、小波分析、偏微分方程分析、线性矩阵不等式控制工具箱等多个模块，并有良好的界面操作功能和非常优秀的使用指导说明。

1．2 天津市强震动台网简介

天津市强震台网是在“九五”实时数据传输强震台的基础上，在“十五”期间重新建设起来的，总共110个强震台（工程力学研究所的SLJ加速度仪暂停用），涵盖了美国凯尼、瑞士GeoSig、SysCom、港震GSMA－2400IP四种类型的仪器，平均间距14．5 km。其中有31个台（光纤通讯）具有准实时波形回传的条件，另外79台当触发条件满足后能在15分钟内自动回传数据到天津市地震局数据处理中心数据库。基于Oracle数据库、Java运行环境及 Web管理的强震动台网实现了强震动数据自动汇集、存储，准实时或定时监控仪器运行状态和台站远程控制与管理等功能。强震台网运行以来，获得了2008年5月汶川M8．0大地震较好的加速度记录和2010年4月9日丰南M4．1地震加速度记录。但此前烈度速报效率低下的情况比较严重。

2 天津市强震动台网烈度速报的几个技术难点及Matlab的应用

通过分析研究，归纳天津市强震动台网速报的技术难点如下：

（1）数据格式不统一，需要厂家的界面化软件对事件文件进行一个一个人工交互式处理（除个别厂家仪器数据能批处理外）。

（2）数据文件内容存在坏死情况，软件无法打开，用解码程序运行则显示错误。

（3）数据文件命名不规范，数据处理后存放路径不统一。

（4）用 Web界面调用Oracle数据，操作不便利，影响下载速度。

（5）数据预处理操作及峰值提取无法批处理。

（6）等值线绘制方法单一且没有融合理论烈度。

以上过程严重影响了烈度速报的速度和效率。针对以上难点，本文利用Matlab处理了下面几个问题：

（1）Matlab对Oracle数据库的直接调用

这个部分主要解决Web调用的数据不能快速下载的缺陷。由于对Oracle数据文件命名作了新的设定（图1），将触发事件和仪器标定事件的命名做了区分，因此可以通过时间段和事件命名的差异直接提取出Oracle某个时间段的地震触发事件文件或标定文件。其涉及的Matlab代码有：

当在本地机器安装好Oracle客户服务端后，配置Matlab对Oracle的驱动，通过以上代码就可以直接从远程Oracle数据库批量读取数据文件。

（2）Matlab与exe文件的融合

该部分主要是充分利用厂家提供的exe格式的解码程序，对从Oracle数据库读出的数据文件进行批量转换成ascii码。由于台站代码和命名固定，所以直接通过文件名就可以确定采用哪个厂家的exe文件进行数据转换。命名方法见图1，如果是标定文件则在上述文件名的基础上加前缀FT＿，即function test的缩写。

Matlab与exe的调用可以通过Matlab提供的system或dos两个命令来操作．由于存在数据文件内容坏死而无法解码的情况，因此可以进一步使用Matlab提供的try．．．catch．．．end的命令来尝试解码，当解码不成功时直接跳过，从而保证后续程序的连续运行。

图1 强震事件文件命名规则Fig．1 Naming rule of strong motion record file．

（3）Matlab的批量显示与通道选择

这部分可通过Matlab界面化函数uicontrol（gcf，‘style’，‘slider’，‘unit’，‘normalized’，‘position’），uicontrol（gcf，‘style’，‘listbox’，‘Position’）组合的方式，提供一个左侧的滑动选择框及右侧的图形曲线显示框（图5）来确认哪道数据可以使用。

（4）仪器烈度的换算

该部分先对第（3）部分选择出来的数据道保存，取每道触发点前20s的数据均值作为基线，以去除各道基线漂移。然后利用ButterWorth带通滤波器压制低频和高频的部分波形能量，考虑强震仪器的宽频与高通特性以及建筑物反应谱设计曲线［10］，本文通频带设为20Hz～2s。之后程序将提取每个台站两个水平道上的第20个最大点作为有效峰值，按每套仪器满量程时对应的gal值，计算出每个仪器两水平道第20个最大点对应的gal值的均值。提取绝对值排序后的第20个最大值作为峰值，可以有效避免仪器尖脉冲的影响（单道记录有时存在多个尖脉冲干扰），同时因为采样率为100sps，不会严重漏用记录峰值及其附近区数据。程序经多次测试，该选点一般可取第5～25点，换算出的gal值较稳定。然后按中国地震烈度表 GB／T 17742－2008［11］中烈度与gal的数值表，采用线性内插的方式（Matlab命令：interp1）获取gal对应的仪器烈度值。当然interp1命令还提供了最近邻、三次多项式、样条法等多种插值方法。由于烈度表不提供V以下的PGA参考值，因此这些外推的低烈度区仪器等值线仅作参考，且其本身因烈度低，其灾评意义亦不明显（Ⅵ以上才有明显地震灾情）。另外，本文亦尝试利用美国ATC－3规范中的加速度有效峰值反应谱EPA（线性化方法计算程序见附录）转换仪器烈度，但由于地震样本量太少，尚未能深入分析本地区PGA，EPA的具体差异。

（5）Matlab的仪器烈度显示

对于烈度的显示，除了利用Delaunay三角剖分三角形显示之外［9］，本部分还提供了两种绘制方案．一是利用Matlab的scatter命令来对不同空间点上的烈度值按数据大小决定显示点大小及颜色的离散方式显示（图7（a）），这种方式比等值线法更能体现单点的仪器烈度。二是利用华北地区理论烈度衰减模型（需要人工输入地震震中的经纬度、震级和长轴方向这些参数，图4（c））构建大网格点（本地区取0．2°）的理论烈度场，然后将仪器记录烈度值加密并取代附近的理论烈度，之后利用griddata命令对非规则的加密后的数据插值，生成烈度等值线。这种方法既利用了仪器烈度也使用了理论烈度，是一种综合效果（图7（b））。而地理底图的显示则直接利用Matlab与arcgis地图格式兼容的命令shaperead和mapshow进行地理图读取与显示。当然，低烈度区（IV以下）的理论计算值亦仅作为画图参考，其本身灾评意义较弱。

（6）可移植性问题

考虑程序的可移植性，程序使用过程中所涉及的程序路径都是相对路径，不存在放在不同路径无法运行的问题。其次，对于台站名命名、地理底图、厂家提供的解码程序等静态资料，程序统一从一个cfg文件夹获取，cfg随整个程序绑定，只要修改cfg的内容并按图1的命名方法对数据文件命名，且替换相应的地理底图，则本文编写的程序可以移植到其他省台网使用。同时在最后烈度成图的操作上本文采用了调用计算过程中生成的文件内容方式，因此当有野外实地烈度调查结果时可以补充到这个烈度结果文件中，从而进一步突出烈度等值线对灾情现场烈度的反映。

利用Matlab人工处理天津市强震台网数据与烈度显示流程及使用的Matlab主要命令如图2所示。

（7）“事件侦听”

由于地震是突发事件，为了保证烈度速报的产出速度，本文的Matlab程序除了提供人机交互处理的方式外，还提供了自动触发处理方式——“事件侦听”（图4）。这是建立在天津市强震动台网数据上传的具体特征的基础上。一般地，地震到来后在5～13分钟内触发的台站基本把数据自动上传到了中心数据库。而无地震的情况下，上传到数据库的文件在5～13分钟内一般小于2个，因此可以每5分钟对数据库上传事件的个数进行统计，如果超过一定的限度，则认为有丛集事件记录上传，并按如下4个步骤进行自动处理（图3）。

图2 Matlab处理天津市强震动台网数据流程图Fig．2 Flow chart of dealing with the data of Tianjin strong motion network using Matlab．

图3 自动扫描处理示意图Fig．3 Sketch of auto－scaning and auto－processing．

①如果有地震事件发生，下次循环时间由5分钟改为4．5分钟并等待8分钟，预留0．5分钟提供程序自动处理前14分钟的数据，并且提取满足扫描条件时刻的前1分钟的数据，因此总共处理8＋1＋5（其中5＝4．5＋0．5）分钟内上传到数据库里的数据。当没有重复事件在13（8＋5）分钟内发生时，数据将提取14分钟的内容，与上一个5分钟的循环有1分钟的重叠。目前我们设定5分钟内上传数大于5个则触发，这个约束条件将可能丢失4个文件，因为1分钟内最多可以上传4个文件。提前一分钟的重叠，至少可以挽回2个文件。

②如果事件发生后，在13分钟之后又发生一次地震事件（即第13～18分钟再次出现丛集事件上传，则提取这个后续事件的前14分钟的事件文件附加到上一次事件中，此时有1．5分钟的重叠）。总共处理14＋14－1．5＝26．5分钟内上传的文件数据。

③如果在18分钟后才再次发生地震事件丛集上传，则作为第二个地震事件处理。

④如果不发生其他事件，则正常进入下一个5分钟的循环。

经过一年的实践运行结果显示，未出现误自动处理现象。当“设定5分钟内上传数大于2个”时，偶尔出现自动处理结果，这些结果常是由仪器工作不正常，或职工当时进行仪器维修造成。

3 应用实例

本实例采用2010年4月9日丰南M4．1地方震加速度记录，触发台站36个，可利用的数据台27个（由于当时外围网络数据传输程序有缺陷，SYSCOM公司的仪器记录文件坏死情况较严重）。地震发生后，当时数据处理过程耗用了多个小时。本程序编写完成后，在3分钟内就可以通过界面化（图4～6）的鼠标操作完成了所有的操作，其中点击“人工数据提取”可以弹出图4（b）的参数输入界面。数据绘制的烈度显示结果（图7）也与当时各区县群众打电话询问的频度很相近，其中宁河县，塘沽区和市中心高层住户有明显震感。需说明的是，图5、图6的纵坐标做了如下处理：先对各台站记录数据归一化以压制异常记录，然后对每台记录的基线递增进行显示以避免曲线堆叠。数据提取按内存变量值为准。

图4 天津强震动台网速报软件界面（a），获取事件选择框（b）和理论烈度计算参数输入框（c）Fig．4 GUI of intensity rapid report of Tianjin strong motion（a），GUI of selecting events（b）and GUI of input parameters to calculate theoretical intensity（c）．

4 认识与讨论

本文将Matlab软件提供的调用数据库命令，调用仪器厂家的exe命令，数据分析处理命令和界面化操作等命令进行有机组合，从而比较成功地解决了天津市强震动台网在烈度速报中遇到的速度与效率问题，大大缩减了仪器烈度发布的时间，由此得到以下认识：

图5 Uicontrol函数实现的数据显示及通道选择操作Fig．5 Data imaging and selecting with uicontrol function．

图6 Uicontrol函数实现的数据显示局部放大（由右滑键控制）Fig．6Local zooming－in with uicontrol function（Operated by right－side slider）．

（1）Matlab简明的程序代码书写格式和丰富的命令式程序库，可以有效降低应用程序开发难度，大大缩减程序代码量。但Matlab软件再强大，也需要和研究者所研究对象的物理模型、数学模型充分融合才能发挥其作用。本文中涉及的理论烈度计算的非线性迭代方程，数据ButterWorth滤波频段设计，“事件侦听”的循环判断方案等等都是具体物理模型的数学化，是Matlab本身不具备的。

图7 丰南4．1级地震的仪器烈度散点显示（a）和等值线显示（b）Fig．7 Equipment intensity imaging by scattering（a）and contouring method（b）for M4．1earthquake at Fennan county．

（2）仪器烈度与传统烈度的“多因子综合判定”、“宏观分级”、“以结果表述原因”三个特征尚有一定的衔接困难，本文的PGA直接转化仪器烈度亦有较多值得商榷之处，利用EPA方法尚还在尝试，由于缺乏足够的地震样本未能深入研究其与PGA的差别。日本目前已不再过分强调仪器烈度与传统烈度的灾害评估差异，并减少现场宏观调查的工作量，突显仪器烈度在快速灾评的指导，他们这种处理方式不失为一种明智之举［3，12］。

由于在本文中并未对天津地区场地进行烈度校正，且没有进一步讨论利用反应谱的结果来刻画仪器烈度的可行性，因此仪器烈度速报的结果还是比较粗糙。另外，Matlab的通讯工具箱直接对串口实时数据，对数字采集器直接访问的研究也未深入展开，在人机交互处理界面上也不太美观，这些工作都是后续需要努力的。

致谢：本文程序编写过程中得到了长安大学邓亚虹博士，北京航空航天大学吴鹏博士的帮助，在此非常感谢！

［1］ 李山有，金星，马强，等．地震预警系统与智能应急控制系统研究［J］．世界地震工程，2004，20（4）：21－26．

［2］ 金星，廖诗荣，陈绯雯．区域数字地震台网实时速报系统研究［J］．地震地磁观测与研究，2007，28（1）：64－72．

［3］ 张敏政．地震烈度及其评定［J］．防灾科技学院学报，2010，12（1）：1－6．

［4］ 金星，张红才，韦永祥．基于地震台网资料快速发布的震动烈度标准及其应用研究［J］．国际地震动态，2008，（10）：20－27．

［5］ 李俊，苏枫，米宏亮，等．ShakeMap及其在地震动快速预估中的应用［J］．中国地震，2010，26（1）：103－111．

［6］ D J Wald，V Quitoriano，T H Heaton，et al．TriNet“ShakeMaps”：Rapid Generation of Peak Ground Motion and Intensity Maps for Earthquakes in Southern California［J］．Earthquake Spectra，1999，15（3）：537－555．

［7］ K T Shabestari，F Yamazaki．A proposal of instrumental seismic intensity scale compatible with MMI evaluated from three－component acceleration records［J］．Earthquake Spectra，2001，17（4）：711－723．

［8］ 泽仁志玛，陈会忠，何加勇，等．震动图快速生成系统研究［J］．地球物理学进展，2006，21（3）：809－813．

［9］ 刘双庆，谢静，聂永安，等．Delaunay三角剖分算法在地震烈度速报中的应用－－以天津市强震台网为例［J］．地震地磁观测与研究，2011，32（5）：115－122．

［10］ 国家标准建筑抗震设计规范管理组 编．建筑抗震设计规范（GB50011－2010）统一培训教材［M］．北京：地震出版社，2010．

［11］ 中国国家标准化管理委员会．中国地震烈度表GB／T 17742－2008［S］．北京：中国标准出版社，2009．

［12］ 颜心仪．利用台湾宽频地震网从事强震预警研究［D］．台湾省：台湾大学地质科学研究所，2007．

附录：Newmark方法计算地震动反应谱子程序（Matlab版）：