MATLAB软件在鹤岗地震台地震观测记录中的应用

第３ｌ卷第４期　防灾减灾学报　Ｖ０ｌ＿３１　Ｎｏ．４　Ｄｅｃ．２０１５　２０１５年ｌ２月　ＪＯＵＲＮＡＬ　ＯＦ　ＤＩＳＡＳＴＥＲ　ＰＲＥＶＥＮＴＩＯＮ　ＡＮＤ　ＲＥＤＵＣＴＩＯＮ　文章编号：１６７４—８５６５（２０１５）０４—００５６—０４　ＭＡＴＬＡＢ软件在鹤岗地震台地震观测　记录中的应用　姜　博，胡宝慧，常金龙，李大伟，张　浩，张　震　（鹤岗地震台，黑龙江鹤岗１５４１０１）　摘要：随着数字化地震计的应用和发展，数字观测资料中不仅记录了宽频带地震信息，也包含了大量的　干扰信号。在实际工作中，对地震信息的分析处理影响很大，有时还会造成地震波形记录不完整或地震　信号淹没在干扰信号中，很难再进行深入的分析和应用。因此，利用ＭＡＴＬＡＢ软件设计ＩＩＲ数字滤波器　对不同记录中的干扰波进行滤波器设计，达到消除干扰波的目的。实例表明，此方法在台站应用效果良　好。　关键词：数字地震记录；干扰波；滤波器；ＭＡＴＬＡＢ　中图分类号：Ｐ３１５．３；Ｐ３１５．８　文献标志码：Ａ　ＤＯＩ：１０．１３６９３／ｊ．ｃｎｋｉ．ｅｎ２１—１５７３．２０１５．０４．０１０　矩阵运算、信号处理和图形显示于一体，工　０　引言　鹤岗地震台经“十五”数字化改造后，　使用甚宽频带ＣＴＳ—ｌ型地震计，极大的丰富　了地震观测信息，同时台站记录的地震信号　常常叠加有某些频率的复杂的背景信号，这　具箱中包含各种经典和现代的数字信号处理　技术，能实现各种数字滤波器的设计口】。文中　给出了ＩＩＲ数字滤波器设计原理和ＭＡＴＬＡＢ　程序设计方法，并给出了用这种方法滤除数　字地震记录中干扰波的应用实例。　些信号就是所谓的噪声［１－２］，这些噪声严重影　响了地震观测质量，给我们的日常分析带来　了很多麻烦。现在几乎所有的台站都使用宽　频带或甚宽频带地震计，这些地震计可以记　１　数字滤波器设计原理　数字滤波器是数字信号处理技术的重要　内容。和模拟滤波器一样，数字滤波器的主　录到各种信号，如地脉动，海浪，爆破，汽　车干扰等。还有通常在较长周期的远震记录　要功能是对数字信号进行处理，最常见的处　理是保留数字信号中的有用频率成分，去除　信号中的无用频率成分　】。针对此特点，我　们对数字地震观测记录中的地震信号和干扰　信号进行分离，并最大限度的保留地震信号，　去除干扰信号。根据鹤岗地震台的观测记录　上叠加有高频振动；在近震记录上叠加有远　震面波或由于天气的变化形成的长周期的信　号［３－４］ｏ　ＭＡＴＬＡＢ是一套用于科学计算的可视　化高性能的语言和软件环境，它集数值分析、　收稿日期：２０１５—０９—２０　修订日期：２０１５—１０—２５　作者简介：姜博（１９８４一），男，本科，助理工程师，现主要从测震研究方面的Ｔ作。Ｅ—ｍａｉｌ：ａ３０５５９５７＠１２６一ｍ　４期　姜博，等：ＭＡＴＬＡＢ软件在鹤岗地震台地震观测记录中的应用　５７　特点，运用ＭＡＴＬＡＢ程序设计了无限脉冲响　应ＩＩＲ数字滤波器，根据滤波器的性能指标　要求，设计滤波器的分子和分母多项式系数。　数字滤波器的设计目的是使滤波器的频率特　性达到所给定的性能指标。其性能指标也包　括通带波纹、阻带衰减、通带边界频率、阻　带边界频率等ｎ】。　ＩＩＲ滤波器的传递函数为：　酢　＝　：鬈㈩　ｈ（ｎ１为滤波器的脉冲响应，ｎ＝Ｏ～∞均有　值。Ｍ和Ｎ为分解的分子和分母多项式的系　数个数。理想数字滤波器的频率特性在通带　内必须满足：　ｊｌ　Ｊ　　（　）＝一ａｏ９　Ｋ　（２）　数字滤波器的设计　ＩＩＲ经典设计就是将已设计好的模拟滤波　器按一定变换原理转换为数字滤波器。ＩＩＲ数　字滤波器经典设计法的一般步骤为：　（１）根据给定的性能指标和方法不同，　首先对设计性能指标中的频率指标，如数字　边界频率进行变换，转换后的模拟频率指标　作为模拟滤波器原型设计的性能指标。　（２）估计模拟滤波器最小阶数和边界频　率。　（３）设计模拟低通滤波器原型。利用　ＭＡＴＬＡＢ工具函数ｂｕｔｔａｐ，ｃｈｅｂｌａｐ，ｃｈｅｂ２ａｐ，　ｅｌｌｉｐａｐ等。　（４）由模拟原型低通滤波器经频率变换　获得模拟滤波器（低通、高通、带通、带阻　等）。　（５）将模拟滤波器离散化获得ＩＩＲ数字滤　波器。　根据以上步骤，下面我们以设计　Ｂｕｔｔｅｒｗｏｒｔｈ滤波器为例，介绍具体的设计步骤　和方法。　（１）确定滤波器的性能指标，主要包括：　采样频率Ｆｓ、通带边界频率Ｆｃｐ、阻带边界　频率Ｆｃｓ、通带波纹Ｒｐ、阻带衰减Ｒｓ。　（２）用函数【Ｎ，Ｗｎ］＝ｂｕｔｔｏｒｄ（ｗｐ，ＷＳ，Ｒｐ，Ｒｓ）　求得数字滤波器的最小阶数Ｎ和归一化截止　频率Ｗｎ。　（３）用函数【ｂ，ａ］＝ｂｕｔｔｅｒ（Ｎ，ｗｎ）求出滤波　器传递函数多项式系数向量ａ和ｂ。　（４）用函数［Ｈ，ｆ］＝ｆｒｅｑｚ（ｂ，ａ，Ｎｎ，Ｆｓ）求得滤　波器的频率特性。　（５）用函数Ｙｔ＝ｆｉｈｅｒ（ｂ，ａ，Ｘｔ）对输入信号进　行滤波。　（６）用命令ｓｕｂｐｌｏｔ（２，１，１），ｐｌｏｔ（ｔ，Ｘｔ）；ｔｉｔｌｅ（’　输入信号’）绘制输入信号图，ｓｕｂｐｌｏｔ（２，１，２），　ｐｌｏｔ（ｔ，Ｙｔ），ｔｉｔｌｅ（’输出信号’）绘制输出信号图。　ＭＡＴＬＡＢ滤波程序的代码如下：　ｌｏａｄ　ｔｅｘｔ．ｔｘｔ；％调入信号数据文件　Ｘｔ＝ｔｅｘｔ；％原始波形数据　Ｆｓ＝５０；％采样频率５０Ｈｚ　ｗｐ＝Ｆｃｐ　２，Ｆｓ；ｗｓ＝Ｆｃｓ￥２／Ｆｓ；％转换为归一　化频率　Ｒｐ＝１；Ｒｓ＝２５；％通带衰减和阻带衰减　Ｎｎ＝１２８；％绘频谱图所用点数　［Ｎ，Ｗｎ］＝ｂｕｔｔｏｒｄ（ｗｐ，ＷＳ，Ｒｐ，Ｒｓ）；％求滤波器　的最小阶数和临界频率　［ｂ，ａ］＝ｂｕｔｔｅｒ（Ｎ，Ｗｎ）；％求Ｂｕｔｔｅｒｗｏｒｔｈ数字　滤波器传递函数的分子ｂ和分母ａ　ｉｆｇｕｒｅ（１）；％图形（一）　［Ｈ，ｆ］＝ｆｒｅｑｚ（ｂ，ａ，Ｎｎ，Ｆｓ）；％滤波器特性图　ｓｕｂｐｌｏｔ（２，１，１），ｐｌｏｔ（ｆ，２０　ｌｏｇｌ０（ａｂｓ（Ｈ）））；　ｘｌａｂｅｌ（’频率／Ｈｚ’）；ｙｌａｂｅｌ（’振幅／ｄＢ’）；ｇｒｉｄ　ｏｎ；　Ｓｕｂｐｌｏｔ（２，１，２），ＰｌＯｔ（ｆ，１　８０／ｐｉ　Ｕｎｗｒａｐ　（ａｎｇｌｅ（Ｈ）））　ｘｌａｂｅｌ（’频率／Ｈｚ’）；ｙｌａｂｅｌ（’相位／＾ｏ’）；ｇｒｉｄ　ｏｎ；　ｎ＝０：ｌｅｎｇｔｈ（Ｘｔ）一１；ｔ＝ｎ／Ｆｓ；％转换为时间序　列　ｉｆｇｕｒｅ（２）；％图形（二）　ｓｕｂｐｌｏｔ（２，１，１），ｐｌｏｔ（ｔ，Ｘｔ）；ｔｉｔｌｅ（’输入信号’）；　％绘制输入信号　ｘｌａｂｅｌ（’时间／ｓ’）；ｙｌａｂｅｌ（’振幅’）；％坐标轴　标识　Ｙｔ＝ｆｉｈｅｒ（ｂ，ａ，ｘｔ）；％对输人信号进行滤波　５８　防灾减灾学报　３１卷　ｓｕｂｐｌｏｔ（２，１，２），ｐｌｏｔ（ｔ，Ｙｔ），ｔｉｔｌｅ（’输出信号’）；％　绘制输出信号　ｘｌａｂｅｌ（’时间／ｓ’）；ｙｌａｂｅｌ（’振幅’）；％坐标轴　标识　３数字滤波器运行实例　随着鹤岗地震台周边的环境的变化，仪　器记录到的波形叠加了很多低频和高频的成　分，对地震记录影响很大，尤其对中小地震　影响很大，很难分析，为了把其中的干扰信　号排除，针对不同的频率成分，设计了不同　时阊，ｓ　鹤岗地震台滤渡后渡彩　图１　原始记录波形和滤波后波形　Ｆｉｇ．１　Ｏｒｉｇｉｎａｌ　ｒｅｃｏｒｄｉｎｇ　ｗａｖｅｆｏｒｎｌ　ａｎｄ　ｉｆｌｔｅｒｅｄ　ｗａｖｅｆｏｒｍ　如图３可以看出鹤岗地震台记录到的地　方震还叠加了许多低频干扰，所以设计通带　边界频率为０．５Ｈｚ，阻带边界频率为０．２Ｈｚ，　通带波纹为ｌｄＢ，阻带衰减为３０ｄＢ（图　鹋韵地煮台ｌ原始记录波彤　善　州　．Ｉ　枷蝴：　矗　一　１　出Ｉ　～　．一　～　一　骧　ｆ　｜１．　…一　　－图３原始记录波形和滤波后波形　Ｆｉｇ．３　Ｏｒｉｇｉｎａｌ　ｒｅｃｏｒｄｉｎｇ　ｗａｖｅｆｏｒｍ　ａｎｄ　ｉｆｌｔｅｒｅｄ　ｗａｖｅｆｏｒｍ　的数字滤波器进行滤波。　鹤岗地震台记录到的波形经常叠加了许　多高频干扰，正常波形受到严重干扰，由图１　可以看出记录波形中有很多高频干扰，因此　可以运用低通滤波器将高频干扰滤除，所以　设计通带边界频率为４Ｈｚ，阻带边界频率为　８Ｈｚ，波纹为ｌｄＢ，阻带衰减为２５　ｄＢ（图２）。　使代码中Ｆｃｐ＝４，Ｆｃｓ＝８运行后得图１，　可以看到低于２Ｈｚ的信号可以通过，高于　８Ｈｚ的信号不能通过，通带波纹为ｌｄＢ，阻带　衰减为２５ｄＢ，由图２可见符合设计要求。　图２频率相应　Ｆｉｇ．２　Ｆｒｅｑｕｅｎｃｙ　Ｒｅｓｐｏｎｓｅ　４），使代码中Ｆｃｐ＝０．５，Ｆｃｓ＝０．２，［ｂ，ａ］＝ｂｕｔｔｅｒ　（Ｎ，ｗｎ，’ｈｉｇｈ’），运行后如图３，小于０．２Ｈｚ的　信号完全被滤除，可以清晰的分析出震相　、　Ｓｇ。　图４频率相应　Ｆｉｇ．４　Ｆｒｅｑｕｅｎｃｙ　Ｒｅｓｐｏｎｓｅ　４期　姜博，等：ＭＡＴＬＡＢ软件在鹤岗地震台地震观测记录中的应用　５９　以上实例都是运行文中的ＭＡＴＬＡＢ代码　最合适的边界频率，尽量完整保留原始地震　实现的，只是根据不同的干扰频率特性，设　信息。　计了不同的通带边界频率Ｆｃｐ和阻带边界频　随着社会的发展，地震台站周围的观测　率Ｆｃｓ，随着不同干扰频率的改变，我们可以　环境发生了巨大的改变，干扰信息越来越多，　设计的Ｆｃｐ和Ｆｃｓ也作相应的改变，即可达到　对于去除干扰信号，保留完整、准确的地震　滤除干扰信号的目的，方法灵活实用。有利　记录信息还有很大难度，希望通过对地震信　于提高地震波分析的能力和准确性，更提高　号和干扰信号不断的进行分析和积累，来增　了地震观测资料的质量，在台站更加实用。　强台站工作人员分析、识别、处理地震波中　干扰信息的能力。　４结论与讨论　参考文献：　随着科学技术的不断发展，地震计所记　【ｌ】万永革．数字信号处理的ＭＡＴＬＡＢ实现【Ｍ］．北京：科　录到的地震波形信息不仅大大增加了，各个　学出版社，２０１２．　频率的干扰信息也随之增加了很多，而数字　［２】陈建军，许玉红，李兴坚，等．数字滤波器在高台地　信号处理技术正好解决了地震波数据处理方　震台地震记录日常分析中的应用『Ｊ］．西北地震学报，　面的许多问题　】。针对每个台站的观测记录有　２００７，２９（４）．　不同的特点，在排除数字地震信号中的干扰　［３】曾庆堂，起卫罗，马志刚，等．ＭＡＴＬＡＢ消除腾冲台　信号时，通过对台站的干扰信号和正常背景　数字地震记录中干扰波的应用［Ｊ】＿华南地震，２０１４，１　噪声的分析，准确选取通带边界频率和阻带　（３４）．　边界频率，应用适合台站自身特点的滤波函　［４］焦旭霞，马克博，韩英王，等．ＭＡＴＬＡＢ在地震记　数进行滤波，就可以将原始观测记录中的干　录日常分析中的应用初探【Ｊ】＿地震地磁观测与研究，　扰信号尽量的滤掉。同时，滤波后的地震信　２０１０，１（３１）．　号会影响地震震相的准确度和震源参数分析　［５】刘瑞丰，陈翔，沈道康，等．宽频带数字地震记录震　的准确性，因此，我们在设计波器时要选择　相分析【Ｍ】．北京：地震出版社，２０１４．　ＡＮＡＬＹＳＩＳ　ＡＮＤ　ＲＥＳＥＡＲＣＨ　ＯＦ　ＳＥＩＳＭＩＣ　ＯＢＳＥＲＶＡＴＩＯＮ　ＤＩＳＴＵＲＢＡＮＣＥ　ＩＮ　ＨＥＧＡＮＧ　ＪＩＡＮＧ　Ｂｏ，ＨＵ　Ｂａｏ－ｈｕｉ，ＣＨＡＮＧ　Ｊｉｎ－ｌｏｎｇ，ＬＩ　Ｄａ－ｗｅｉ，ＺＨＡＮＧ　Ｈａｏ，ＺＨＡＮＧ　Ｚｈｅｎ　（Ｈｅｇａｎｇ　Ｓｅｉｓｍｉｃ　Ｓｔａｔｉｏｎ，Ｈｅｉｌｏｎｇｊ　ｉａｎｇ　Ｈｅｇａｎｇ　１　５４　１　０　１，Ｃｈｉｎａ）　Ａｂｓｔｒａｃｔ：Ｗｉｔｈ　ｔｈｅ　ａｐｐｌｉｃａｔｉｏｎ　ｏｆ　ｄｉｇｉｔａｌ　ｓｅｉｓｍｏｇｒａｐｈｓ，ｏｂｓｅｒｖａｔｉｏｎｓ　ｎｏｔ　ｏｎｌｙ　ｒｅｃｏｒｄｓ　ｔｈｅ　ｂａｎｄ　ｂｒｏａｄｂａｎｄ　ｓｅｉｓｍｉｃ　ｉｎｆｏｒｍａｔｉｏｎ，ｂｕｔ　ａｌｓｏ　ｃｏｎｔａｉｎｓ　ａ　ｌａｒｇｅ　ｎｕｍｂｅｒ　ｏｆ　ｉｎｔｅｒｆｅｒｉｎｇ　ｓｉｇｎａｌｓ．Ｉｎ　ｐｒａｃｔｉｃｅ，ｔｈｅ　ａｎａｌｙｓｉｓ　ａｎｄ　ｐｒｏｃｅｓｓｉｎｇ　ｏｆ　ｉｎｆｏｒｍａｔｉｏｎ　ｏｎ　ｓｅｉｓｍｉｃ　ｇｒｅａｔ，ｓｏｍｅｔｉｍｅｓ　ｒｅｓｕｌｔｉｎｇ　ｉｎ　ｓｅｉｓｍｉｃ　ｗａｖｅｆｏｒｍ　ｒｅｃｏｒｄｓ　ａｒｅ　ｉｎｃｏｍｐｌｅｔｅ　ｏｒ　ｓｕｂｍｅｒｇｅｄ　ｓｅｉｓｍｉｃ　ｓｉｇｎａｌ　ｉｎｔｅｒｆｅｒｅｎｃｅ　ｓｉｇｎａｌ，ｖｅｒｙ　ｄｉｆｆｉｃｕｌｔ　ｔｏ　ｃａｒｒｙ　ｏｕｔ　ｉｎ－ｄｅｐｔｈ　ａｎａｌｙｓｉｓ　ａｎｄ　ａｐｐｌｉｃａｔｉｏｎｓ．Ｔｈｅｒｅｆｏｒｅ，ｔｈｅ　ｕｓｅ　ｏｆ　ＭＡＴＬＡＢ　ｓｏｆｔｗａｒｅ　ｄｅｓｉｇｎ　ＩＩＲ　ｄｉｇｉｔａｌ　ｆｉｌｔｅｒ　ｆｏｒ　ｄｉｆｆｅｒｅｎｔ　ｒｅｃｏｒｄｉｎｇ　ｉｎｔｅｒｆｅｒｅｎｃｅ　ｗａｖｅ　ｆｉｌｔｅｒ　ｄｅｓｉｇｎ，ｔｈｅ　ｐｕｒｐｏｓｅ　ｏｆ　ｅｌｉｍｉｎａｔｉｎｇ　ｉｎｔｅｒｆｅｒｅｎｃｅ　ｗａｖｅｓ．Ｅｘａｍｐｌｅｓ　ｓｈｏｗ　ｔｈａｔ　ｔｈｉｓ　ｍｅｔｈｏｄ　ｗｏｒｋｓ　ｗｅｌ１．　Ｋｅｙ　ｗｏｒｄｓ＂ｄｉｇｉｔａｌ　ｓｅｉｓｍｉｃ　ｒｅｃｏｒｄｉｎｇ；ｉｎｔｅｒｆｅｒｅｎｃｅ　ｗａｖｅ；ｆｉｌｔｅｒ；ＭＡＴＬＡＢ