地震模型1.doc

地震模型摘要地震预报是当代自然科学领域里一个直接关系到人民生命安全和社会经济的发展，是科学难度很大的前沿课题。近几十年来，人民的努力虽然取得了丰硕的成果，预报的实践有些有限的成功。在地震发生前有很多前兆性指标，如磁场强度、氡值、电压、雨量、水温等，这些指标都从不同的侧面反映了地震活动的各项特征。在正常情况下，也常常有些指标有明显的异常，而令一些指标并不出现异常。这些都给实际的前兆指标数据特征分析以及地震预报工作带来了困难。本文针对地震影响因素多，数据多，联系弱的特点，分别建立了针对各个任务的数学模型，首先，处理初始数据时选取日平均变化来消除一天中数据的随机因素的影响，然后更具这些数据建立了地震发生前后各个指标的不同程度度量模型解决了地震对指标的敏感度的度量问题，进而找出了EW，气氡，水位，电压，雨量这些衡量地震的主要因素，并分析出了这些指标的重要程度，在解决找统计量时，利用上次任务所得的主要因素的组合来衡量一次地震的地震前兆指标数量特征，并得出了电压达到26.8754，气氡达到17.685左右时就有地震发生的可能性，考虑全面周到，效果较好，在模型改进中所提出的判别分析法，科学有效，对数据利用率较高。而除此之外的其他几个指标与震级有关，这样就这些指标的分析在任务三中分别建立了地磁前兆异常的动态从属函数模型以及非线性综合模型，进而找出了地震级数与这些因素之间的非线性关系，对于上述数学方法在地震预测预报中的应用没有达到预期效果，须进一步研究改进。综合这些工作，我们应用易语言编写了分析地震数据的平台，通过分析地震数据生成了报告。 一、 问题重述地震是地壳快速释放能量过程中造成的振动。虽然预测地震是世界性难题，但迄今科学界普遍认为，有可能反映地震前兆特征的指标可能不少于10个。已经有专业仪器在多个定点实时按秒记录这些指标的数据，期望通过对记录数据的分析研究找到地震的前兆特征。现已采集到某地2005年1月1日至2010年6月30日按小时观测的10多个指标的数据，和该地区该时期内已发生地震的时刻、经纬度、震级及震源深度的数据。这些数据中隐藏着地震发生的前兆特征。科学地截取这些数据的有用片段，对数据进行合理地预处理，用数学方法揭示地震前兆的数据特征，是一项很有意义的研究工作。题给数据中的这10多个指标，究竟哪些与地震的发生有关，有何种关系，是单一关系还是复合关系；除这10多个指标外还有哪些因素及含题给指标在内的哪些指标的哪种数学模型更能反映地震的前兆特征等等，人们迄今仍不很清楚，需要进行深入地研究。地震数据的观测是持续进行的，随着时间的推移数据的规模会不断扩大。从中挖掘地震的前兆特征，必须有合理的数学模型，也必须有科学高效的算法分析平台。因此，请参赛队结合附件中给出的实际记录数据，尝试完成以下任务。任务一：分析数据特征，建立数学模型以度量各指标对地震发生的敏感程度。任务二：构造由某些或全部指标构成的综合指标，使其尽可能地集中反映地震发生前的数据特征的统计规律。任务三：结合题给数据，广泛查阅与地震相关的其它指标的数据和分析方法，建立数学模型来研究地震发生前的数量特征。任务四：将前述各项任务的计算程序集结成地震数据分析平台，使其能够完成形如题给数据的其它地震数据的分析，并能自动输出前述任务的重要的分析结果。任务五：对于进一步的研究设想写一篇切实可行的报告。二、问题分析地震是地下岩层受应力作用错动破裂造成的地面震动，是一种破坏性极强的自然灾害，是自然灾害之首恶。地震可以在很短的时间内使一座城市夷为平地，使无数的家庭支离破碎。为了有效地减小地震带来的损失，对地震进行预测则很重要。本题中给出了地震发生前2005年1月1日至2010年6月30日按小时观测的10多个指标的数据和该地区该时期内已发生地震的时刻、经纬度、震级及震源深度的数据。这些数据中隐藏着地震发生的前兆特征。在地震发生前有很多前兆性指标，如磁场强度、氡值、电压、雨量、水温等，这些指标都从不同的侧面反映了地震活动的各项特征。在正常情况下，也常常有些指标有明显的异常，而令一些指标并不出现异常。这些都给实际的前兆指标数据特征分析以及地震预报工作带来了困难。21、问题的性质本文主要是对所给数据进行处理，用数学方法揭示地震前兆的数据特征。2.2、解决问题的难点和关键1、如何确定前兆个指标对地震发生的敏感程度。2、如何构造地震发生的前兆综合指标，尽可能的揭示地震发生前的数据特征的统计规律。3、构建地震数据分析平台。三、模型假设1、假设地震发生前2005年1月1日至2010年6月30日按小时观测的10多指标的数据是准确无误的。2、假设由于其他因素而引起10多个指标数据的变化以及非正常波动可以忽略不计。3、假设第二问剔除由第一问求出的敏感度较低的指标对地震的敏感度可以忽略不计。4、假设地震的前兆指标的数据特征符合一定的概率统计分布。四、符号说明第次地震的地震前的检测指标随时间的变化规律 第次地震的地震后的检测指标随时间的变化规律 第次地震的地震后的预测指标随时间的变化规律地震级别组成的集合地震前兆测量指标集合各因素所占的权重地震前兆预测异常特征量观测值的自相关系数观测曲线对应于时间轴的斜率五、模型的建立与求解5.1数据的处理与分析5.1.1、各指标数据的分析 在地震活动性分析中有许多前兆性的数据指标，例如题中数据给出的氡值、水位、磁场强度等因素。这些指标在不同的侧面反映了地震活动的各种特征。但在实际的预报中，常常有些参数在一些中强以上的地震前出现比较明显。在正常情况下，也常常有些参数出现较明显的异常，而另一些参数并不出现异常。这些都给实际预报带来了困难。5.1.2数据的预处理由于题目中给出地震每小时测数据受其他因素受天气、气候等其他外在因素的影响造成的数据波动较大，为了剔除偶然因素的影响，本文采取了求日平均的方法来对数据进行预处理。如下图分别列举了2005/2006的EW人均值图像。由于篇幅有限，其他年份指标的图像都附在附录里面。图1： 2005年EW日平均变化趋势图图2： 2006年EW日平均变化趋势图本数据预处理使用了均值结构模型消除了其他偶然因素的影响，使地震前兆指标数据更加准确的反应地震发生的前兆信息。5.2模型的建立与求解5.2.1指标敏感度模型设变量的变化规律如下所示： （1）求前后与的不同敏感程度。现已知一函数序列，表示第年的的变化规律，现在结合来求，现在若，现在若已知，如何求出与的关联程度。若关联程度较大，则与不同程度表示：即定义与的不同程度为： （2）分析变量与的关联程度 （3）则由分析可得该问题便为如何求的问题。现已知与，如何根据已知的以上信息求出，已知的变化规律，可知令 （4）模型求解：由上述结论我们分析数据得出各个因素在地震发生前后的敏感程度分别如下：电压EWNS地温水位气温气压水温气氮雨量0.21.80.80.12595.32.34.340.22.80.5表1：各因素与地震敏感度的表但是考虑地震因素影响上述指标外，可能还有其它的因素影响上述指标。其中查阅资料可得出：水位，水温，气压，气温受到外界影响最大，这样我们除上述几个指标剩下的指标基于与地震关联度从大到小排列：气氮EWNS雨量电压地温2.81.80.80.50.20.1 表2：主要地震前兆指标敏感度的表标准化后有：气氮EWNS雨量电压地温0.450.290.130.080.030.01表3：主要地震前兆指标敏感度归一化表如下图分别画了敏感度较高的EW与敏感度较低的低温年平均曲线，在途中可以看出在前兆指标因素发生剧烈变化时，会有地震的发生。图3： 2005，2006年EW日平均变化趋势图图4： 2005、2006年地温日平均变化趋势图5.2.2综合指标模型方法及步骤用两个集合，一是地震级别组成的集合，用表示，二是地震前兆车辆指标集合，用来表示，且有，每个指标因素都有个状态级，如此需要对个指标影响程度分敏感、一般、不敏感、很不敏感五个状态。有个指标因素分别用来表示，其中，于是可得到敏感指标分析数学模式 （5）有了这个模型，分别给予地震级别的隶属度函数值再根据指标因素间的关联度及其重要性，分配权数这样就得到了集合D上的模糊关系： （6）基于地震的综合指标步骤有下列步骤：（1）、详细分析影响地震发生的因素，并筛选出若干重要因子作为分析指标。（2）、在分析地震各指标以及波形的基础上，综合得到地震的等级指标。（3）、给出适当的隶属度和权数。（4）、修正矩阵U。确定分析指标地震发生前的检测指标较多，而各个指标的敏感程度不同，经过第一问的各因素的关联度分析，筛选出气氮、EW、NS等六个指标。隶属度和权数的确定应用综合指标进行判别，隶属度函数是一个关键，用它来反映地震发生。设是评定集合上的一个模糊子集，则隶属度函数为。为隶属度，其值可以有专家评判方法来确定。假设等级划分与指标间基本呈线性关系，则根据第一位求出的关联度来确定隶属度，因此选用，为了方便起见将隶属度扩大一百倍后可得：，权数的确定是人给定的，选取的方法有很多，本文选用特费尔法，并确定气氮、EW、NS、雨量、电压地温的权数分配为：这里为上的模糊子集，按模糊数学的常规技法 （7）模糊关系矩阵的确定 根据，就可以得到在集合上所需的模糊关系矩阵： （8）由第一问可得各因素的敏感度为：电压EWNS地温水位气温气压水温气氮雨量0.21.80.80.12595.32.34.340.22.80.5 表4：地震前兆指标敏感度的表每次地震的等级：D=2.3 3.6 3.1 3.1 3.2 3.0 4.4;稳定性分析数学模式：0.07615.08500.32900.09203.23490.12620.05891.34221.27210.05482.21430.55030.43200.90640.56010.16433.75110.98950.35760.87870.60750.11053.25580.23480.04561.85951.67410.12241.80060.27300.07150.81960.40190.01722.73850.89490.07190.68060.12460.00483.09550.0307表5：每次地震前兆稳定性系数地震指标权重 -6.1024 电压 - 0.5740 东西磁场强度 1.0214 南北磁场强度 2.0882 水位 1.6259 气氡 -0.8900 雨量假设地震发生等级与地震发生前数据关系为D=Ub+g;RI=7337;4948;1932;2259;3265;4473;7401;在7次地震中，选中的每个指标在地震发生前的平均值：26.86113.77374.566717.748812.0109026.9050-2.03324.733617.58668.88190.027727.0236-3.02343.337917.568012.0243027.0298-2.76233.613217.571212.4378026.76480.13136.972117.618833.2737026.7152-1.36346.478917.620217.96150.660226.70211.87578.584417.655027.21010表6：每次地震前每个指标在地震发生前的平均值在这7次中上述每个指标的均值为：26.8574 -0.4859 5.4695 17.6241 17.6857 0.0983在这7次中上述每个指标的方差为：0.1369 2.5458 1.9290 0.0631 9.1568 0.2480由此得出第2,9个指标与地震的级数有相当大的关系，而其他的指标为稳定指标。即当电压达到26.8754，气氡达到17.685左右时就有地震发生的可能性。而1,3,4,10这几个指标与震级有关;假设是线性关系。即：D=a1*E(1)+a2*E(3)+a3*E(4)+a4*E(10),有上述地震发生前的平均值可得系数为：a1 a2 a3 a4=3.1119 0.4654 -4.6999 -0.4878;5.2.3地震前兆指标数量特征模型 先查阅与地震相关的其他指标如震级，频率等，而后先分析单个指标的数量特征，并以电磁波EW为例，利用动态从属函数建立模型，对于气温，水位等其他指标，类推可得。而后利用震级，频率等指标的数学研究方法，将已知条件的10个指标综合起来，利用非线性组组合建立模型。电磁波EW的数量特征研究。地震活动是一种复杂现象，各种观测数据之间缺乏对比性，难以利用统一的物理场进行描述。为了将各种观测量的异常信息进行表达，可把各前兆异常量转化为无量纲量以便进行对比。各种地震前兆观测量虽然是不同的物理量，但其共同的特点都是随时间变化的数值，即都是时间函数，各种地震前兆异常也都表现为各种观测量随时间的突出变化，异常形态虽多种多样，但究其本质共同点都是观测曲线随时间的斜率变化。因此，观测曲线随时间的斜率变化将是判断异常的重要指标。由以上分析从各种前兆数据中提取地震信息量的方法应具以下特点：（1）、地震信息是一个无量纲量；（2）、用前兆观测量Y（t）对时间轴的斜率变化K，反映观测值的速率变化，一般情况下，K值越大，异常量越大；（3）、用前兆观测值的自相关系数r来显示观测值质量的好坏，r值越大，观测数据质量越高。原理和方法 从地震单项前兆观测值中识别地震异常，并计算它属于地震异常的程度，需要遵循的基本条件是：（1）、观测本身内在质量要好；（2）、“异常”出现的频次要低，且出现“异常”与“地震”之间对应程度要高；（3）、“异常”与“正常”相比，有明显的差异性，且界限要清楚。按照地震前兆异常的普适性和上述基本条件，提出如下地震前兆异常特征量： （9）式中： （10） （11）是前兆观测随时间变化的数列，是相应的时间序列 公式（1）中的k物理意义是观测曲线Y(t)对于时间轴的斜率，反映观测速率的变化；r是观测值的自相关系数，反应观测值内在质量的好坏。从（1）式可知，只有k和r都大时，才可能大。在内取值，表示在所研究的资料在时段属于地震异常的程度。作为模糊界点，只有当时，才视为异常，反之，即为正常状态。是一个经过学习得到的常数。它的选择原则和方法如下：（1） 使得“异常”频次降低；（2） 使“异常”和“地震”之间由最好的相关性；（3） 在满足上述两个原则下要尽可能的小；（4） 选择的“地震”样本要视监视区的活动水平和预报目标而定；（5） 选择后的值，应保持稳定，只有出现预报目标与实际情况有明显差异时，才按照新的样本要求进行修改，以保持第二条原则的实现。模型的检验选择样本从属函数中的的选择由样本学习而决定，不同学科不同地域选择的参数可能会不同，需要在实践中不断的完善，一些前兆场可直接根据原始观测数据分析其变化，而地磁场本身是一个变化场，且磁层、电离层强烈扰动时引起的地磁场变化非常剧烈，远大于地震异常引起的变化。从属函数法以曲线斜率作为判断异常的组要指标，如果直接分析某一点的从属函数值，会频繁出现“异常”，判断食物。同样是前兆手段，地磁学在应用从属函数方法时首先对原始地磁场资料预处理，才由利于寻找震磁关系。电磁波由地磁场产生，而地磁场值包含核远源场，外源场及地下局部磁异常源场，提取电磁波EW的问题是如何找出因地磁场的地下局部磁异常发生变化的问题，并以此判断该变化量是否是由地震引起的局部异常。核源场是一个相对稳定场，外源场的变化受诸多因素影响，而且其变化幅度远大于因地震引起的异常量，因此应用在应用资料预报之前，首先应剔除外源场引起的地磁场变化。任一点的地磁场垂直强度Z由以下三部分组成： （12）C-地核源场，e-外源场，N-地下局部磁场。由此可推得下面公式： （13）下面采用均值法消除偶然因素如太阳黑子等造成的突变，并用对应一年中的均和几年的均值作差以消除外源场和核源场的变化而引起地磁场的突变 （14）由此可推得： （15）正常情况下地点确定后局部磁场在很小的范围内波动；由于地震孕育而导致局部磁场发生变化时，就在孕震区产生局部磁异常，这个范围比较小，即使两地都处于孕震区，异常将由于距离震源距离不同而有所差别，最终导致的变化， 代表局部磁异常的相对变化。从属函数的计算 对应年份观测数据和各年份数据的平均作差消除基本场与外空场而引起的EW变化,用于消除外源场和核源场的变化而引起地磁场的突变，取各年份的均值是为偶然因素而引起的电磁波突变经过多年样本学习，确定值，研究电磁波随时间的变化情况。每次地震每个指标相对应的值电压EWNS地温气氡雨量第一次地震0.50000.98670.96660.99940.99880.9916第二次地震0.98830.50000.96690.50000.99970.8144第三次地震0.99840.96130.92550.71230.99370.9236第四次地震0.97860.99160.98820.98890.99990.8648第五次地震0.99940.98370.50000.99231.00000.9836第六次地震0.99270.99500.99430.96130.99991.0000第七次地震0.96300.99890.99490.99820.50000.5000 表7：每次地震前兆预测异常特征量值的分析：的取值范围为，0.5为判别地震发生的界点。当时，则地震没有发生的可能。当时，地震有发生的可能，值越大，则地震发生的可能性越大。每个指标的值 1.0e-003 * 0.0117 0.1812 0.3049 0.0021 0.0142 0.0736图5： EW两年相对变化趋势图5.2.3地震前兆指标数量特征模型针对地震预测研究中，变量的指标较多，数据携带的信息较弱且预测变量是非线性关系，因此造成预测结果稳定性不强的问题，构造了非线性组合预测模型，解决地震预测的上述问题，应用于本题中的地震前兆指标数量特征。即使数据携带的信号很弱，但是我认为这些数据足以通过统计方法来探求深层次的东西，因为数据不服从标准的平稳时间序列和线性模型，所以我们要小心谨慎地选择统计方法。下面构造的多变量非线性模型在一定程度上解决了上述问题且在地震预报应用中取得了较满意的效果。非线性组合模型非线性组合模型是将多个看用于预测的一元线性模型进行线性组合，组成非线性组合模型，其优点是可以增加预测稳定性。非线性组合模型的应用背景设是预测变量，是影响的指标。（1）、指标数据携带信号较少、较弱，需要选择携带信号较强的指标用于预测。 （2）、对于选定的指标需要建立较好的非线性预测模型。 （3）、为增强预测的稳定性，最后构建非线性组合模型用于预测。最佳指标和最佳曲线选择对于每一个指标，对应预测变量，预先设计若干非线性模型，并进行曲线拟合，根据统计检验和拟合优度选择最优指标和最佳曲线。设所选的最优指标为相应的最佳曲线函数是，是最优指标和最佳拟合曲线个数。非线性组合模型 （16）其中是预测变量，是选择出的最佳非线性预测曲线，是自变量指标，是相应的权重，。非线性组合模型参数的估计一般参数估计的方法可选择：（1）等权平均法； （2）加权最小二乘法；（3）方差-协方差方法；（4）最小二乘法；非线性组合模型的应用在本题中，选择最优指标和最佳曲线，通过最小二乘法估计最佳曲线的参数，并构建非线性组合模型。最佳指标和最佳曲线的选择 初步选择影响最大的若干地震参数指标等，见表一为选择最优指标和最佳拟合曲线设计了七种曲线模型，见表二。选择拟合方式以及统计检验，选出6个最优指标以及他们相应的最佳曲线。见表三。2.2、非线性组合模型 考虑到增加预测的稳定性，参数估计采取等权平均法，经过Matlab编程得出了地震各因素的权重如下：电压EWNS地温气氡雨量-0.0000-0.00290.00010.0057-0.0014-2.3620 画出部分地震前兆指标与地震发生的关系图如下：图6：相对震级图5.2.4地震数据分析平台的设计把前三问的程序集成，基于易语言编程，系统开发使用Windows2000server、MS-SQLServer2000数据库，对地震目录、前兆观测数据、分析报告等地震信息进行了系统管理，组成了地震数据分析平台，本平台在实现了地震数据的数字化和自动化管理功能的基础上，系统运用数据挖掘的方法对地震预报进行了尝试，使其能够完成形如题给数据的其它地震数据的分析，并能自动输出前述任务的重要的分析结果。其主界面如下图.图7：地震数据分析平台主界面基于此地震数据分析平台，分析了2009的地震指标数据，得出的结果是地震发生时间为10月31日，震级为4.1级，在误差范围内，实现了较好的地震前兆数据分析与得出重要的分析结果。附录1：日平均值的程序 %每年数据的平均值RI=7344 4968 1944 2280 3288 4488 7416;year=2005 2006 2007 2007 2008 2008 2009;S = load(matlab.mat);er=zeros(40*7,6);for i=1:7 u=strcat(S.t,num2str(year(i),(,num2str(RI(i)-40*24+1),:,num2str(RI(i),1 2 3 4 9 10); data=eval(u); temp1=zeros(40,6); for j=1:40 temp1(j,:)=mean(data(24*j-23):(24*j),:); end er(40*i-39):(40*i),:)=temp1;endK1=zeros(7,6);%每次地震每个指标的K值r1=zeros(7,6);%每次地震每个指标的r值for i=1:7 K1(i,:)=sum(1:40)*sum(er(40*i-39):(40*i),:)-40*sum(diag(1:40)*er(40*i-39):(40*i),:); K1(i,:)=K1(i,:)/(sum(1:40)2-40*sum(1:40).2); r1(i,:)=sum(diag(1:40)*er(40*i-39):(40*i),:)-(1/40)*sum(1:40)*sum(er(40*i-39):(40*i),:); r1(i,:)=r1(i,:)./(sqrt(sum(1:40).2)-(1/40)*sum(1:40)2)*sqrt(sum(er(40*i-39):(40*i),:).2)-(1/40)*sum(er(40*i-39):(40*i),:).2);endgh=zeros(7,6);for j=1:6 gh(:,j)=K1(:,j).*r1(:,j);endmiu=zeros(7,6);%每次地震每个指标的miu值for j=1:7 miu(j,:)=1./abs(1+min(gh)./gh(j,:);end2年平均曲线的Matlab程序function s=years_avg()yinsu=dataVIEWNSETWHATAPWTAIDRAINPNSPEW;year=2005 2006 2007 2008 2009 2010;qs=input(请输入你所需要的第几年的编号：);qw=input(请输入您要的因素编号：);coona=database(sql server,sa,xuenhappy);setdbprefs(DataReturnFormat,numeric);s=zeros(1,372);for i=1:12 for j=1:31 if i10&j10 t=strcat(yearqs,-0,num2str(i),-0,num2str(j),%); t=strcat(char(39), t,char(39); t= ,t; else if i9 t=strcat(yearqs,-0,num2str(i),-,num2str(j),%); t=strcat(char(39), t,char(39); t= ,t; else if i9&j10 t=strcat(yearqs,-,num2str(i),-0,num2str(j),%); t=strcat(char(39), t,char(39); t= ,t; else t=strcat(yearqs,-,num2str(i),-,num2str(j),%); t=strcat(char(39), t,char(39); t= ,t; end end end s(j+31*(i-1)=fetch(coona, strcat(use earthquck select avg(,yinsuqw,) from years,yearqs,$ where Convert(Varchar(10),data,120) like,t) ); if isempty(s(j+31*(i-1) s(j+31*(i-1)=0; end ende