基于matlab的地震活动性分析

Matlab在地震活动性图像分析中的应用1），李红光2）1）河北省地震局2）中国地震应急搜救中心摘要：地震活动性分析是地震预测、地震工程的一个重要依据，地震活动性分析又多是通过图像来表现。

Matlab是一种简单易学、强大的计算功能和编程可视化的计算机语言。

本文用Matlab语言编程，实现了地震统计区内地震的快速选取，并根据这些选中的地震进行地震活动性分析。

关键词：Matlab语言；地震活动性引言地震活动性研究就是通过分析一定震级区间内的地震时间、空间的分布特征，探讨其物理含义，进而对地震发生的规律进行科学总结。

通过地震活动性研究，可对地壳介质非均匀性和运动形态有宏观的了解和总体把握，因此可服务于地震预测和地球动力过程等研究[1]。

在地震安全性评价中，通过地震活动性分析，为工程场地一定时间内的地震活动性趋势和地震环境做出评价，为合理划分潜源区和确定其地震活动性参数提供依据[2]。

地震活动图像的分析方法很多，有简单的图像描述法，如地震震中分布、蠕变曲线、M-T图等；也有采用统计参数表征地震活动时空图像特征的方法，如b值、地震活动度S、地震能流密度、地震强度因子MF分布等。

Matlab具有强大的计算能力、计算结果可视化和编程效率高的优势，它是地震活动性分析的一个有力工具。

Matlab是1984年由美国MathWorks公司推出的荣誉产品。

早在20世纪80年代中期，Matlab就在我国出现，大规模流行时再90年代中期以后。

现在Matlab已被广泛应用在科学研究、工程计算等方面。

M atlab采用全新的数据类型和面向对象编程技术，采用了新控制流和函数结构，特别是包含很多常用的子函数，非计算机专业人员非常容易用Matlab来实现很复杂的计算程序。

并且Maltab提供了图像处理功能，可以很方便的生成图形。

在地震活动性数字图像分析中，用Matlab可以很简单、方便的实现研究人员的思想。

Matlab在地震预测、地震工程、地震波记录处理和地震定位等方面得到了应用。

王凤基于Matlab的BP预测模型进行地震前兆预测的研究[3]；焦旭霞利用Matlab分析地震日常记录[4]；李晋平利用Matlab 实现了地震数据的小波变换[5]；李敬利用Malab中的数字信号处理工具箱进行数字地震记录中干扰波的去除[6]。

Matlab以它的简单易学、强大的工具箱和结果可视化等优势已经在地震行业被广泛应用。

本文通过编写Matlab程序，快速画出任意地震统计区M-T图和应变释放曲线图来说明Matlab在地震活动性分析中的优势。

一计算地震与地震统计区关系画M-T图和应变释放曲线图首先要选择地震统计区内的所有符合要求的地震目录。

地震统计区是CPSHA方法中特有的概念，区内大小地震的发生具有统计上的协调性，满足G-R关系。

而且地震统计区的边界是决定地震样本的不同归属，具有很重要的意义[7]。

所以在地震安全性评价中为了准确的统计处地震统计区地震活动性参数，我们必须根据地震统计区的边界精确的判定出地震的归属。

我们可以把地震统计区看作一个任意的多边形，地震震中作为一个点，判断地震是否发生在地震统计区内，可以转换为判断一个点是否在多边形内的问题。

用计算机实现这个问题，不仅大大的加快了工作效率，还提高了结果的精度。

二判断点与多边形的关系理论点与多边形位置关系的判断问题是计算机图形学中的基本问题，也是多年来国内外的诸多学者比较关注的问题，他们针对点在多边形内的检测问题作了广泛的研究，提出了许多算法以提高计算效率[8]。

归纳起来常用的方法有3种：叉积判断法、夹角之和检验法与射线法．其中，叉积判断法要同时对全部边进行叉积运算，角度法则要使用复杂的三角运算，因此二者的计算量都非常大．目前，工程上应用最多的还是射线法，这种方法简单可靠。

本文采用射线法来解决这个问题，下面简单的介绍射线法的原理。

地震统计区是一组首尾相连的线段的有限集合，线段的端点就是地震统计区的顶点，所以用有序的顶点序列P１，P２，⋯ Pｎ，Ｐｎ＋１表示地震统计区，其中Ｐｎ＋１＝Ｐ１.每两个相邻的点连接起来就是地震统计区的一个边。

射线法判断点与多边形关系的关键是通过判断点引一条射线，求射线与多边形边的交点个数，如果交点个数为偶数则点在多边形外，如果为奇数则点在多边形内．如图1所示，有5个判断点分别对应5条水平向右的射线，它们与多边形相交的位置各不相同：如①-⑥只与多边形的边相交；②-⑦和③-⑧与多边形的顶点相交；射线④-⑨和⑤-⑩与多边形的某些边重合．不同情况下交点个数的求解方法不同，下面分析不同情况下交点个数的求解过程．１.一般情况如图1中的①-⑥只与多边形的边相交，不通过多边形的顶点，属于射线与边相交的一般情况，此时可以很容易得到交点个数。

２．射线通过多边形顶点的情况。

射线通过多边形顶点的情况有2种：I、通过单个或多个互不相邻的顶点；II、通过2个或多个相邻顶点，即射线与多边形的边重合。

如图１中所示的射线②-⑦和③-⑧属于第1种情况，即射线②-⑦通过单个顶点，射线③-⑧通过2个不相邻顶点。

射线②-⑦通过多边形的顶点P２，该顶点对应的2条边都位于顶点的上侧；射线③-⑧通过多边形的顶点P１１和P４其中P１１对应的2条边都位于顶点的下侧，P４对应的2条边则分别位于顶点的两侧．由图1可以看出：射线②-⑦通过顶点P２并与边P１P１２和P３P４各有一个交点，如果把P２算作一个交点，则射线与多边形的交点个数等于3，根据射线法的奇偶判断原则，点②位于多边形内。

但从图１中可以看出点②位于多边形外，因此判断有误；射线③-⑧通过P４和P１１顶点，点③应该在多边形外，而从图1可看出，点③在多边形内，由此看出，当射线通过多边形的顶点时，交点个数的计算与射线通过的顶点所对应边的情况有关。

当射线通过多边形的顶点时，如果该顶点的前后2个相邻顶点位于射线的同侧，则认为在相交的顶点处有2个交点或者没有交点，因为断定点是否在多边形内的主要依据是射线与多边形边交点个数的奇偶性，所以认为这种顶点不作为交点，如顶点P２。

如果该顶点的前后2个相邻顶点位于射线的异侧，则认为在相交处的顶点就是一个顶点，如顶点P４。

射线通过相邻顶点的情况如图１所示的射线④-⑨和⑤-⑩，这2条射线分别通过多边形边P５P６和P８P９，此时由于线段与射线重合，所以可以分别把P５、P６和P８、P９各作为一个交点，此后的交点个数的判断与射线通过单个顶点时的判断方法相同。

3.判断点在多边形的边上的情况。

如果判断点与边的下侧端点连线的斜率等于该边斜率，则判断点在多边形的边上。

此外，如果判断点与多边形的某个顶点重合，判断点在多边形的边上。

图１任意多边形与任意点的几种关系根据以上理论，我们给出图2判断点与多边形的关系流程图：三利用Matlab绘制M-T图和应变释放曲线与分析根据以上判断点与多边形关系的方法编制程序，并根据台网中心提供的地震目录和第四代区划图提供的地震统计区划分，我们可以快速的绘制任意地震统计区内的地震活动性分析图像。

我们以Ｍ－Ｔ图和应变释放曲线图为例子，来表明本方法可以为工程地震提供很大的方便。

我们选取的华北平原地震统计区作为例子。

图3是根据本方法得到的华北平原地震统计区的Ｍ－Ｔ图和应变释放曲线图，由于华北地区(除黄海及边远地区外)Ms≥4¾地震自1484年之后基本完整[9]，所以我们取起算震级为Ms=4¾，起始时间为1484年。

图3 华北平原地震统计区M-T图和应变释放曲线图从图3中可以看出：自1485年以来华北平原地震带有两个地震活跃期（1485－1730年和1791年开始的活跃期），每个活动周期按地震累积和应变释放又可分为两个活动段。

第一活动周期的第一活动段从1485年河北遵化5级地震到1624年滦县6½级地震，属于前兆释放阶段，最大地震是6½级。

第二活动段是1658年涞水6级地震，一直到1730年“北京北郊”6½级地震，是属于大释放阶段，最大地震是1679年三河—平谷8级地震。

第二活动周期的第一活动阶段是1820年许昌6级一直到1882年深县6级地震，是属于加速释放阶段。

最大震级是磁县7½级地震。

第二活动段是1937年菏泽7.0级地震一直到1983年菏泽6.0级地震是属于大释放阶段，最大地震是唐山7.8级。

这个活动周期内还发生1830年磁县7½级地震。

目前可能处于该地震带的第二活动周期的后期，今后百年内该带属于地震剩余应变释放阶段或下一个活动期的应变积累阶段。

估计未来地震活动水平低于其长期平均水平。

四总结与讨论：（1）本文用Matlab语言编程，根据点与多边形的关系判断理论，实现了地震统计区内地震的快速选取，并根据这些选中的地震进行地震活动性分析。

（2） Matlab是一种简单易学、具有强大的计算功能和结果可视化的计算机语言。

它被广泛应用在工程计算等方面，近年来也被应用到地震数据分析和地震工程计算等领域。

（3）根据作者自己编写的程序，可以非常方便和快捷的画出任意地震带任意时间段的地震M-T图和应变释放曲线图。

Matlab绘制的图像不仅是矢量的，而且图中的其他属性都可以方便任意改变。

（4）本次工作只是利用Matlab对地震活动性进行了简单的分析，以Matlab同样可以应用的地震分析的其他领域。

参考文献：[1] 王建.地震活动性研究及其应用于地震预测的一些问题[J].中国地震，2005，21（3）：451～456[2] 胡聿贤. 地震安全性评价教程[M].1999.北京：地震出版社[3] 王凤等. 基于MATLAB的BP预测模型在地震前兆预测中的应用研究[J]. 华北地震科学，2009，27（1）：48～51[4] 焦旭霞等. MATLAB在地震记录日常分析中的应用初探[J]. 地震地磁观测与研究，2010，31（1）：118～122[5] 李晋平等. 用MATLAB实现地震数据的小波变换[J]. 化探物探计算技术，2002，24（2）：163～168[6] 李敬等. 数字地震记录中干扰波的排除[J]. 防灾技术高等专科学校学报，2004，6（3）：20～25[7] 潘华，金严等地震带与地震统计区关系探究［Ｊ］. 地震学报，2003，125（3）：308～313[8] Feito F R，Torres J C，Urena L A．Orientations i mplicity and inclusion test for planar polygons[J]．Computer& Graphic，1995，19(4)：595~600[9] 黄玮琼，李文香等，中国大陆地震资料完整性研究之一——以华北地区为例［Ｊ］，地震学报，1994，16(3)：273～280附录：绘制地震带的M-T图和应变释放曲线图Matlab程序代码%m-tclear all;mm=input('请输入起算震级:');Byear=input('请输入起始时间:');if (mm>10)|(mm<1.0)fprintf(1,'输入的震级不正确\n');return;endif (Byear>2050)|(Byear<-1000)fprintf(1,'输入的时间不正确\n');return;end%Byear是统计地震起始时间；%读入地震带数据；'地震带序号：\n'cd E:\m-ttu\mid格式;evd=dir('*.mif');lev=length(evd);for ii=1:levllname(ii)=length(evd(ii).name)-4;ddname(ii,1:llname(ii))=evd(ii).name(1:(strfind(evd(ii).name,'.')-1)); sss=strcat(int2str(ii),',', ddname(ii,1:llname(ii)),'\n');fprintf(1,sss);enddddd=input('请选择要计算的地震带序号：');if (dddd>lev)|(dddd<1)fprintf(1,'输入的地震带序号不正确\n');return;end% dddd=str2num(ss);% for dddd=1:levfid=fopen(evd(dddd).name,'r');for jj=1:9tline = fgetl(fid);endN1=fscanf(fid,'%f\n',[1 1]);N(dddd)=N1(1,1);xx=fscanf(fid,'%f %f\n',[2 N(dddd)]);xx=xx';fclose(fid);cd E:\m-ttu;filename1='EQ16全国地震目录Ms5.EQT';filename2='mttu.txt';fid=fopen(filename1,'r');fid1=fopen(filename2,'w');nnn=0;while 1% a=fread(fid,'%s %s',[2 1]);tline = fgetl(fid);if ~ischar(tline), break, end%disp(tline)nnl=length(tline);nm=strfind(tline,' ');nt=strfind(tline,'.');% dizhi=native2unicode(tline((nm(2)+1):nnl));ms=str2double(tline((nt(3)-1):(nt(3)+2)));if ms>=mmsitex(1)=str2double(tline((nt(2)-3):(nt(2)+3)));sitey(1)=str2double(tline((nt(1)-3):(nt(1)+3)));%sitex(1)为震中经度；sitey(1)为震中纬度%计算点与多边形的关系:nn=0;for ii=1:(N(dddd)-1)pd=(xx(ii,2)-sitey(1))*(xx(ii+1,2)-sitey(1));%计算两条直线相交情况if pd<0% nn=nn+1;xxb=(xx(ii+1,1)-xx(ii,1))*(sitey(1)-xx(ii,2))/(xx(ii+1,2)-xx(ii,2))+xx(ii,1);if xxb<=sitex(1)nn=nn+1;xx1(nn)=xxb;endendif pd==0if (xx(ii,2)~=sitey(1))&(xx(ii+1,2)==sitey(1))%与多边形上的顶点相交情况xxb=xx(ii,1);if xxb<=sitex(1)if (ii~=N(dddd)-1)&((xx(ii+2,2)-sitey(1))*(xx(ii+1,2)-sitey(1))<0)nn=nn+1;xx1(nn)=xxb;endif (ii==N(dddd)-1)nn=nn+1;xx1(nn)=xxb;endendelseif (xx(ii,2)==sitey(1))&(xx(ii+1,2)==sitey(1))nn=nn+1;%水平射线与多边形边重合endendendif mod(nn,2)==1% disp(tline(1:nm(2)))% disp(dizhi)% nnn=nnn+1;year=str2double(tline(2:5));month=str2double(tline(6:7));ri=str2double(tline(8:9));fprintf(fid1,'%d %d %d %f %f %f\n',year,month,ri,sitey(1),sitex(1),ms); endendendfclose(fid1);fclose(fid);% end%画图fid=fopen('mttu.txt','r');[a,count]=fscanf(fid,'%d %d %d %f %f %f\n',[6 inf]);a=a';%大概的换算为年，然后吧年作为横坐标a(:,1)=a(:,1)+a(:,2)/12+a(:,3)/365;fclose(fid);% plot(a(:,1),a(:,6));for ii=1:(count/6)if (a(ii,1)>=Byear)btt=ii;break;endendb=a(btt:(count/6),6);b(1)=sqrt(10^(4.8+1.5*b(1)));for ii=2:(count/6-btt+1)b(ii)=b(ii-1)+sqrt(10^(4.8+1.5*b(ii)));endb=b/(3.16*10^8);hold on;subplot(2,1,1);gcf1=stem(a(btt:(count/6),1),a(btt:(count/6),6),'k');set (gcf1,'Marker','none');ylim([4 10]);xlim([Byear 2050]);mtname=strcat(ddname(dddd,1:llname(dddd)),'地震带M-T图');title(mtname);ylabel('震级M');xlabel('时间T');subplot(2,1,2);gcf2=stairs(a(btt:(count/6),1),b,'k');text('Interpreter','latex','String','$$\sqrt{E}(3.16\times10^8J^\frac{1}{2})$$','Po sition',[Byear b(count/6-btt+1)],'FontSize',16);ybname=strcat(ddname(dddd,1:llname(dddd)),'地震带应变释放曲线');title(ybname);xlabel('时间T');xlim([Byear 2050]);。