基于概率统计模型的地震数据分析与处理.doc

摘要本文主要以我国华东地区地震活跃带为例，分析了基于概率统计模型的分析方法在地震数据处理中的应用.为了使模型更加科学和简单，我们采用马尔科夫更新过程来模拟地震发生的点过程.经过对地震数据的分析，提出用指数分布与威布尔分布的混合分布来模拟较大地震发生的时间间隔所服从的概率分布，采用极大似然估计方法求出未知参数，并验证了模型的有效性.关键词：马尔科夫更新过程； 指数分布； 威布尔分布； 地震数据处理.Abstract (请校正)In this paper, a seismically active belt in the east China area as an example, analyzes the model based on probability and statistics analysis method in the application of the seismic data processing. In order to make the model more scientific and simple, we use markov renewal process to simulate the earthquake point process. Through the analysis of seismic data, using the mixed distribution of exponential distribution and weibull distribution to simulate the large earthquake obedience by interval probability distribution, using maximum likelihood estimation method to find the unknown parameters, and verify the validity of the model.Keywords: markov renewal process; index distribution; Weibull distribution ;of seismic data processing.i目录摘要iAbstract (请校正)ii第一章 引言1第二章 预备知识22.1马尔科夫更新过程简介22.1.1马尔科夫更新过程22.1.2马尔科夫更新过程推论22.2威尔分布简介3第三章 基于概率统计模型的分析方法在地震数据处理中的应用举例33.1数据处理33.2实验53.2.1震级频率柱状图53.2.2时间间隔的箱线图63.2.3采用标准的指数分布来拟合数据集73.2.4采用标准的威布尔分布拟合数据集83.3地震发生的概率计算和模型验证103.3.1 地震发生的概率计算公式103.3.2 概率计算及分析验证10第四章 总结与展望11参考文献12附录12致谢18iii第一章 引言马尔可夫过程（Markov process）是一类随机过程。它的原始模型马尔可夫链，由俄国数学家A.A.马尔可夫于1907年提出。该过程具有如下特性：在已知目前状态 （现在）的条件下，它未来的演变 （将来）不依赖于它以往的演变 ( 过去 ) 。在现实世界中，有很多过程都是马尔可夫过程，如液体中微粒所作的布朗运动、传染病受感染的人数、车站的候车人数等，都可视为马尔可夫过程。关于该过程的研究，1931年A.H.柯尔莫哥洛夫在概率论的解析方法一文中首先将微分方程等分析的方法用于这类过程，奠定了马尔可夫过程的理论基础.威布尔分布：在可靠性工程中被广泛应用，尤其适用于机电类产品的磨损累计失效的分布形式。由于它可以利用概率值很容易地推断出它的分布参数，被广泛应用与各种寿命试验的数据处理.在概率论和统计学中，指数分布（Exponential distribution）是一种连续概率分布。指数分布可以用来表示独立随机事件发生的时间间隔，比如旅客进机场的时间间隔、中文维基百科新条目出现的时间间隔等等。指数分布可以看作当威布尔分布中的形状系数等于1的特殊分布，指数分布的失效率是与时间t无关的常数，所以分布函数简单。本文着重以地震的数据处理为例探讨基于概率统计模型的分析方法的相关应用.在近几十年，世界上很多国家都建立了基于本土文化的地震危险性评估（SHA）系统.就可考证文献记载，描述地震发生的点过程模型基本分为两个大类：泊松过程和非泊松过程5.在以往的研究基础上，本文选取马尔科夫更新过程67.这种更新过程能够模拟较大地震发生的随机性和复杂性，尽管其本身并不复杂.实际上，我们是希望利用它来捕捉传统的泊松模型所不能捕捉到的信息，同时又使得分析过程简单化.单一使用指数分布或者威布尔分布其效果都不是太好,文章中采用指数分布与威布尔分布的混合分布来提高模型的准确度.第二章 预备知识2.1马尔科夫更新过程简介2.1.1马尔科夫更新过程在有限的状态空间中，表示在正实数上的分布函数，表示在上的转换矩阵，表示在上的概率分布.让我们来考虑一种定义在完备概率空间上的二维随机过程：(1);(2)对于每一个，有;(3)对于每一个且，都有那么就称为由决定的马尔可夫更新过程.2.1.2马尔科夫更新过程推论(1) 是带有转换矩阵和初始分布的马尔可夫链.(2) 当，且时，有条件独立.在本文中，马尔可夫链表示连续可见的状态，过程表示连续的等待时间.在我们的应用中，可视化的状态是已经发生过的地震，可以按照不同的震级来进行分类.2.2指数分布简介定义：若随机变量X的概率密度为：其中（0)为常数，则称随机变量X服从参数为的指数分布。指数分布的累计分布函数为指数分布的特性：指数函数的一个重要特征是无记忆性（Memoryless Property，又称遗失记忆性）。这表示如果一个随机变量呈指数分布当s,t0时有P(Ts+t|Tt)=P(Ts)即，如果T是某一元件的寿命，已知元件使用了t小时，它总共使用至少s+t小时的条件概率，与从开始使用时算起它使用至少s小时的概率相等。2.2威尔分布简介在可靠性工程中被广泛应用，尤其适用于机电类产品的磨损累计失效的分布形式。由于它可以利用概率纸很容易地推断出它的分布参数，被广泛应用与各种寿命试验的数据处理。双参数威尔分布是一种单峰的正偏态分布函数，其概率密度函数表达式为：式中：k形状参数，无因次变量； c尺度参数，其量纲与速度相同。第三章 基于概率统计模型的分析方法在地震数据处理中的应用举例表1：华东地区5级以上地震（1970年-2011年）3.1数据处理在这里我们选取地震局的关于华东地区地震数据，从中得到该地区地震的时间及震级数据.本文主要针对5级以上地震进行研究，从中经过筛选得到关于5级以上地震数据如表1所示：时间（年月日）震级时间（年月日）震级时间（年月日）震级时间（年月日）震级197008105.0197903025.4198605235.1200204225.4197106055.2197906195.5198702175.4200207235.1197112275.0197907096.3198912255.3200209035.0197112305.3197907115.2199002105.5200303235.1197210125.2198003095.0199101295.5200303305.1197308045.1198008025.1199111055.1200311255.0197311305.0198111096.1199201235.6200511266.0197312315.6198202145.8199407265.6200511265.2197404225.8198204225.0199509205.6200607045.5197509025.7198311076.2199611096.4201006055.0197509025.1198405215.6199707285.5201010245.0197610065.8198405216.2199802105.0201101125.2197611025.0198511305.6199809075.0201207205.3根据以上表1所列的数据，可以把相邻两次地震的时间间隔求出来，并列于下表：表2：表1所列地震事件的时间间隔（单位：年）0.82144596651 0.29622336377 0.73443302892 0.25470700152 0.56008751903 0.05670091324 2.85437785388 0.10971080670 0.00845319635 0.00523592085 0.12540525114 0.55511035008 0.78275494673 0.65994101979 0.96923706240 0.02083713851 0.81148021309 0.39717656012 0.76875761035 0.65234018265 0.32077815830 1.26904680365 0.21471461187 2.00218607306 0.08775304414 0.26605022831 2.50776255708 0.00046803653 0.30747526636 010.60628044140 1.36320585997 1.54192732116 13.92043759513 0.00036339422 0.54046423135 0.71650494673 0.38493150685 1.09265601218 0.00000190259 0.53603500761 0.21625000000 0.07196347032 1.52454337900 0.57496385084 1.52315829528 2.33446537291 0.47862823440 3.62268264840 对表2的数据进行排序可得表3：表3：排序之后的时间间隔00.128767120.400000000.715068491.2712329000.183561640.476712330.739726031.3643836000.216438360.536986300.767123291.528767100.005479450.219178080.539726030.786301371.545205500.008219180.252054790.550684930.810958902.005479500.019178080.265753420.561643840.819178082.328767100.054794520.298630140.572602740.967123292.506849300.07397260.306849320.602739731.095890402.854794500.084931510.323287670.657534251.139726003.624657500.115068490.386301370.663013701.153424703.923287703.2实验3.2.1震级频率柱状图根据表1提取的震数据，利用matlab绘制震级频率柱状图如下图所示：（程序见附录，程序1）图1 震级频率柱状图由图1可明显看出，发生5.0级地震次数较多，发生6.0及以上级地震次数较少，发生5.9级地震次数为0即断开点为5.9级，且以此为断开点后震级呈现两段走的趋势.根据对马尔科夫更新过程的理解我们可将前后两个状态分别对应马尔科夫更新过程的两个状态，亦即：将震级,定义为状态“1”，将定义为状态“2”.由此我们可以将华东地区5级以上地震数据分成两部分进行研究.3.2.2时间间隔的箱线图根据表2所提供的数据，我们可以利用matlab绘制时间间隔的箱线图，如下图所示：（程序见附录，程序2，3）图2 时间间隔的箱线图图3 时间间隔的直方图从图2中可以清楚地看出，在所研究的数据集合中存在一些异常值，一般地处理方式是删除这些异常值，但在我们的研究中，从地震孕育的物理背景考虑，我们推测这些异常值的出现告诉我们在所研究的数据集合中可能存在两种不同的行为模式，也就是说存在采用两种分布的混合对数据集合进行拟合的可能性.同时从图3中我们发现地震间隔时间的范围大体在03.9之间，大多数时间间隔小于2.0，部分数据大于2，仅有2个时间间隔值大于3.这样的数据特征反映出在华东地区5级以上地震的时间间隔中表现出双重的行为特征，其中具有较大时间间隔我们采用Weibull分布，对于具有较小时间间隔的数据，通过大量实验，我们采用指数分布来描述，这一点我们可以通过下面3.2.3和3.2.4来说明：3.2.3采用标准的指数分布来拟合数据集我们采用标准的指数分布来拟合表3时间间隔的数据，其结果如下图所示：（程序见附录，程序4）图4：标准的指数分布拟合由图4可知对于较小的时间间隔用指数分布拟合效果较好，而对于较大的时间间隔其效果不是很好.3.2.4采用标准的威布尔分布拟合数据集我们再利用matlab工具采用标准的威布尔分布对表3的时间间隔进行拟合，其结果如下图所示：（程序见附录，程序5）图5：标准的威布尔分布拟合由图5可看出，对于较大的时间间隔用威布尔分布的拟合效果较好.综合3.2.3的结果可以采用一个指数分布与威布尔分布的混合分布来拟合数据.其混合分布的概率密度函数的形式为： （5）其中为指数分布的概率密度函数，为威布尔分布的概率密度函数.形式如下: （6） （7）利用极大似然估计对混合分布的概率密度函数中的参数进行估计，如下表所示：（程序见附录，程序6）表4:混合分布中的参数估计值中的参数中的参数权重=0.8225a=1.5054b=6.5579q=0.9026将上述估计值代入（5）式子，得混合分布的概率密度函数，即我们所研究的对象服从指数分布和威布尔分布的混合分布.由（5）式计算对应的累积分布函数作图如下：（程序见附录，程序7） 图6 实验数据和由（5）式推导出来的理论数据的累计分布图从图6中看出，实验数据的累计分布图和混合分布的累积分布曲线具有很好地一致性，说明华东地区5级以上的地震间隔时间较好地服从了指数分布和威布尔分布的混合分布.3.3地震发生的概率计算和模型验证前面我们由已知的数据，通过分析，建立了数据模型.我们要利用这个数据模型对华东地区未来5级以上地震的发生趋势做出判断，下面我们分为几个步骤来阐述这个问题. 3.3.1 地震发生的概率计算公式地震发生概率的计算问题主要是讨论在已知上一次发生5级以上地震的前提下，经过时间的平静期之后，发生下一次5级以上地震的可能性，即概率值.基于上述考虑，我们给出地震发生的概率表达式： （7）和分别表示最后一次发生的地震和下一次将发生的地震，表示最后一次发生的地震和下一次将发生的地震之间的时间间隔，表示预测的时间跨度.假设为对应的分布函数，同时我们研究的是5级以上的地震，因此对应的马尔科夫更新过程中只有一个状态（），因此（7）式也可表示为： （8）3.3.2 概率计算及分析验证我们以华东地区1970年1月1日-2011年1月31日之间的5级以上地震目录为训练集合，建立服从指数分布和威布尔分布的混分分布函数，由（8）式计算继上一次5级以上地震之后，经过时间的平静期之后（），未来（=1年，2年，3年，4年，5年）之后再次发生5级以上地震的概率，计算结果如表5所示：（程序见附录，程序8）表5：利用（1970年1月1日-2011年1月31日）数据计算华东地区未来发生5级以上地震的概率1年0.64140.77830.70412年0.92050.93440.91233年0.97650.98060.97404年0.99300.99420.99235年0.99790.99830.9977运用我们建立的概率模型，训练样本集中最后一次发生的5级以上地震（2011.1.12，），根据地震目录，我们已经知道的下一次在华东地区发生5级以上的地震是在2012.7.20，.根据我们给出的模型，大致取1.5年，计算出发震概率为0.8514.同时，从表5中的第1列，我们也能看出从开始，从最后一次地震（2011.1.12，）起，未来1年内和2年内发生的概率分别是0.6414和0.9205，未来1年半左右发生5级以上地震的概率接近0.85.上述分析说明我们所采用概率模型在分析华东地区5级以上地震是有一定效果的.第四章 总结与展望本文提出了一种新的模型来对华东地区过去几十年地震发生的时间间隔进行分析，并且对实验结果与过往数据进行了对比验证，发现了该方法与实际数据的吻合度较好.这种方法诉诸马尔科夫更新过程与极大似然估计原理，其本质是一种统计学方法.这里开发的方法所使用的数据没有任何的地球物理方面的资料，这使得这种方法并不能真实的触及地震发生的本质，从而并不能用于精确的地震预测工作.所以，在未来的研究中，我们期待着一种基于本文所提供的方法，同时考虑到地球物理科学的假说来提供精确的地震预测，并且希望得到事实的验证.由于地震专业知识和本人能力水平有限，以及时间有限，本文的研究工作还能进一本深入和提高，本人对此将一些不足和有待完善的研究工作展望如下：（1）时间仓促，只是实现了该模型的粗略预测，并没有得到非常精确的预测模型。（2）对于地震的震级，我们只是选取了华东地区的5级及以上地震作为研究对象，并没有将研究的对象进一步延伸。（3）可以拓宽思路，结合地质学和其他新兴科学中寻找思路，以期能找到很好的地震预测模型，能真正为地震预测服务。相信随着研究的深入和研究人员的不懈努力，最终能找到一种能够服务于社会的地震预测模型，能够为人类安全做出贡献。参考文献1 MATLAB R2007基础教程/刘慧颖 编著.北京：清华大学出版社，2008.72 应用随机过程/张波，商豪编著. 2版. 北京：中国人民出版社，20093 附录程序1： x=5.05.25.05.35.25.15.05.65.85.75.15.85.05.45.56.35.25.05.16.15.85.06.25.66.25.65.15.45.35.55.55.15.65.65.66.45.55.05.05.45.15.05.15.15.06.05.25.55.05.05.25.3;a,b=hist(x,15);a=a/length(x);bar(b,a);ylabel(频率);xlabel(震级)程序2： format longx=0.000001902590.0003633940.0004680370.0052359210.0084531960.0208371390.0567009130.071963470.0877530440.1097108070.1254052510.1869006850.2147146120.216250.2547070020.2660502280.2962233640.3074752660.3207781580.3849315070.397176560.4786282340.5360350080.5404642310.555110350.5600875190.5749638510.6062804410.6523401830.659941020.7165049470.7344330290.768757610.7827549470.8114802130.8214459670.9692370621.0926560121.1377511421.1513184931.2690468041.363205861.5231582951.5245433791.5419273212.0021860732.3344653732.5077625572.8543778543.6226826483.920437595;boxplot(x);ylabel(相邻两次地震发生的时间间隔（单位：年）);xlabel(华东地区)程序3： test=load(C:UsersLiJieDesktopdata.txt);x=sort(test);hist(x,20);程序4： format long;x=0.000001902590.0003633940.0004680370.0052359210.0084531960.0208371390.0567009130.071963470.0877530440.1097108070.1254052510.1869006850.2147146120.216250.2547070020.2660502280.2962233640.3074752660.3207781580.3849315070.397176560.4786282340.5360350080.5404642310.555110350.5600875190.5749638510.6062804410.6523401830.659941020.7165049470.7344330290.768757610.7827549470.8114802130.8214459670.9692370621.0926560121.1377511421.1513184931.2690468041.363205861.5231582951.5245433791.5419273212.0021860732.3344653732.5077625572.8543778543.6226826483.920437595;x=sort(x);probplot(exponential,x)程序5： format long;x=0.000001902590.0003633940.0004680370.0052359210.0084531960.0208371390.0567009130.071963470.0877530440.1097108070.1254052510.1869006850.2147146120.216250.2547070020.2660502280.2962233640.3074752660.3207781580.3849315070.397176560.4786282340.5360350080.5404642310.555110350.5600875190.5749638510.6062804410.6523401830.659941020.7165049470.7344330290.768757610.7827549470.8114802130.8214459670.9692370621.0926560121.1377511421.1513184931.2690468041.363205861.5231582951.5245433791.5419273212.0021860732.3344653732.5077625572.8543778543.6226826483.920437595;x=sort(x);weibplot(x)程序6：test=load(C:UsersLiJieDesktopdata.txt);x=sort(test);c1 = x(4:46,:);c2 = x(47:50,:);c3 = c1;c2;% p1,pa1=mle(c2,distribution,exponential);p1,pa1=mle(c1,distribution,exponential);phat= mle(c2,distribution,wbl);mixedpdf=(x,mu,s1,s2,rho)(rho*exppdf(x,mu)+(1-rho)*weibpdf(x,s2,s2);%p,pci1=mle(c3,pdf,mixedpdf,start,p1,phat(1),phat(2),0.95);p=mle(c3,pdf,mixedpdf,start,p1,phat(1),phat(2),0.6)mu= p(1)s1=p(2)s2=p(3)rho=p(4)x=0:0.05:4y=(rho*exppdf(x,mu)+(1-rho)*weibpdf(x,s2,s