关于地震问题的数学模型探究

关于地震问题的数学模型探究吴瑞刘思谢磊摘要地震定位是地震学中最经典、最基本的问题之一，对于研究诸如地震活动构造、地球内部结构、震源的几何构造等此类地震学中的基本问题有重要意义。此外，基于快速准确的地震定位的地震速报，对于震后的减灾、救灾工作也是至关重要的。因此，地震学家一直在不断改进或提出新的定位方法。地震定位问题的提法如下:根据台站对地震到时的观测资料，来确定震源的空间坐标和发震时刻，有时还给出对解的评价。早期的地震定位方法以几何作图法为主。近三十年来由于计算机技术的飞速发展和广泛应用，基于科学计算和计算机技术的智能化数值自动定位方法也得到了迅速发展，并业己成为当前地震定位的主流方法。我国最初的地震定位工作由李善邦先生于1930年在北京鹜峰地震台开创，1953年开始采用多台站大规模观测数据确定震中，现在大多使用国际流行的定位方法。针对本文的数据处理方面，主要运用MATLAB软件进行数据处理。有关的图像处理，也是运用MATLAB软件进行处理即根据局部地形、地貌等高线图，绘制相应的三维地形、地貌曲面图。关键字:最小二乘法Geiger方法差值拟合P波S波前言1.问题的背景中国位于世界两大地震带一环太平洋地震带与欧亚地震带的交汇部位，受太平洋板块、印度板块和菲律宾海板块的挤压，地震断裂带十分发育。大地构造位置决定，地震频繁震灾严重。在20世纪里，全球共发生3次8.5级以上的强烈地震，其中两次发生在我国;全球发生两次导致20万人死亡的强烈地震也都发生在我国，一次是1920年宁夏海原地震，造成23万多人死亡;一次是1976年河北唐山地震，造成24万多人死亡。这两次地震死亡人数之多，在全世界也是绝无仅有的。上世纪以来，中国共发生6级以上地震近800次，遍布除贵州、浙江两省和香港、澳门特别行政区以外所有的省、自治区、直辖市。中国地震活动频度高、强度大、震源浅，分布广，是一个震灾严重的国家。1900年以来，中国死于地震的人数达55万之多，占全球地震死亡人数的53%;1949年以来，100多次破坏性地震袭击了22个省(自治区、直辖市)，其中涉及东部地区14个省份，造成27万余人丧生，占全国各类灾害死亡人数的54%，地震成灾面积达30多万平方公里，房屋倒塌达700万间。地震及其他自然灾害的严重性构成中国的基本国情之一。统计数字表明，中国的陆地面积占全球陆地面积的十五分之一，即百分之六左右;中国的人口占全球人口的五分之一左右，即百分之二十左右，都不到百分之二十，然而中国的陆地地震竟占全球陆地地震的三分之一，即百分之三十三左右，而造成地震死亡的人数竟达到全球的1/2以上。当然这也有特殊原因，一是中国的人口密、人口多;中国的经济落后，房屋不坚固，容易倒塌，容易坏;第三与中国的地震活动强烈频繁有密切关系。统计，20世纪以来，中国因地震造成死亡的人数，占国内所有自然灾害包括洪水、山火、泥石流、滑坡等总人数的54%，超过1/2。从人员的死亡来看，地震是群害之首;而在经济上所造成的损失，最大的主要是气象灾害(洪涝)，气象我国的地震活动主要分布在五个地区的23条地震带上。这五个地区是:①台湾省及其附近海域;②西南地区，主要是西藏、四川西部和云南中西部;③西北地区，主要在甘肃河西走廊、青海、宁夏、天山南北麓;④华北地区，主要在太行山两侧、汾渭河谷、阴山一燕山一带、山东中部和渤海湾;⑤东南沿海的广东、福建等地。我国的台湾省位于环太平洋地震带上，西藏、新疆、云南、四川、青海等省区位于喜马拉雅一地中海地震带上，其他省区处于相关的地震带上。灾害所造成的经济损失要比地震大的多。图1我国地震带分布图图2中国强震及地震分布图图3中国地震和火山图4中国地震震中分布图2.概念阐述2.1最小二乘法最小二乘法(又称最小平方法)是一种数学优化技术，它通过最小化误差的平方和找到一组数据的最佳函数匹配。最小二乘法是用最简的方法求得一些绝对不可知的真值，而令误差平方之和为最小。最小二乘法通常用于曲线拟合。很多其他的优化问题也可通过最小化能量或最大化嫡用最小二乘形式表达。人们往往对由某一变量七或多个变量f......f构成的相关变量Y感兴趣。如弹簧的形变与所用的力相关，一个企业的盈利与其营业额，投资收益和原始资本有关。为了得到这些变量同Y之间的关系，便用不相关变量去构建Y。尝试用如下函数模型(1)最小化问题的精度，依赖于所选择的函数模型。q个相关变量或p个附加的相关变量去拟和。通常人们将一个可能的、对不相关变量七的构成都无困难的函数类型充作函数模型(如抛物线函数或指数函数)。参数x是为了使所选择的函数模型同观测值Y相匹配。(如在测量弹簧形变时，必须将所用的力与弹簧的膨胀系数联系起来)。其目标是合适地选择参数，使函数模型最好的拟合观测值。一般情况下，观测值远多于所选择的参数。其次的问题是怎样判断不同拟合的质量。高斯和勒让德的方法是，假设测量误差的平均值为0。令每一个测量误差对应一个变量并与其它测量误差不相关(随机无关)。人们假设，在测量误差中绝对不含系统误差，它们应该是纯偶然误差，围绕真值波动。除此之外，测量误差符合正态分布，这保证了偏差值在最后的结果Y上忽略不计。确定拟合的标准应该被重视，并小心选择，较大误差的测量值应被赋予较小的权。并建立如下规则:被选择的参数，应该使算出的函数曲线与观测值之差的平方和最小。用函数表示为:(2)用欧几里德规则表达为:（3）最小化问题的精度，依赖于所选择的函数模型。2.2P波在测震学中，震中距在1000公里以上的称为远震。如果把测量点作为地球的一极，穿过地心和测量点相对应的为另一极，中间的距离按照和地心与测量点之间的张角划分为180度，也可以说震中距在10度以内的称为近震，10度一一105度的为远震，105度以上的称为极远震。对于远震和极远震，不能再象近震一样认为地震波是直线传播的，必须考虑地震波在传播时的折射。从震源发出，一路不断折射传播过来的纵波叫做P波。它是远震测量中的首波。P波(P-waveorprimarywave)是二种体波(体波的命名是因为此波穿越地球内部，相对于体波的是表面波，另一种体波是S波(secondarywave)中的一种。P波意指(primarywave)或是压力波(pressurewave)。在所有地震波中，P波拥有最快的传递度速，因此地震发生时，P波是最早抵达测站，并被地震仪纪录下来的地震波，这也是P波名称的由来。P波的P也能代表压力(pressure)，来自于其震动传递类似声波，属于纵波的一种(或疏密波)，传递时介质的震动方向与震波能量的传播方向平行。对于地球内部构造的了解和推论，大部分是藉由观测地震波中的体波。地震波在不同介质有不同传播时间和路径，在介质交界面时会产生反射、折射，以及相位的改变，地震学家利用这些特性来获得地球内部资讯。当体波穿越地球液态层时，P波在经过下部地函与外地核时会稍许折射。造成P波在104。至1400间会有阴影区，即地震仪记录不到。2.3S波S波(S-waveorsecondarywave)是二种体波(体波的命名是因为波穿越地球内部，相对于体波的是表面波)中之一。它是因地震而产生的，被地震仪记录下来。命名为S波(secondarywave)是因为它的速度仅次于P波(最快的球内部粒子的震动方向与震波能量传递方向是垂直的。s波与P波不同的是，S波无法穿越外地核。所以S波的阴影区正对著地震的震源。2.4走时地震走时是指地震波从震源传到观测点所经过的时间。地震纵波和横波传播的速度不同，同一距离传播时间也不同。传播时间的长短，还受到所经之处的物质组成情况的影响，准确地测出走时，对研究地震有重要意义。3.模型的建立与求解3.1问题一的模型的建立与求解3.1.1问题的重述假定地面是一个平面，在这个平面上建立坐标系见图5。图5中给出了10个地震观测站点(A-J)的坐标位置。图5地震观测站点示意图2009年4月1口某时在某一地点发生了一次地震，图5中10个地震观测站点均接收到了地震波，观测数据见表1。地震观测站横坐标x（千米）纵坐标y（千米）接收地震波时间A50033004月1日9时21分9秒B3002004月1日9时19分29秒C80016004月1日9时14分51秒D140022004月1日9时13分17秒E17007004月1日9时11分46秒F230028004月1日9时14分47秒G250019004月1日9时10分14秒H29009004月1日9时11分46秒I320031004月1日9时17分57秒J34001004月1日9时16分49秒表1地震观测站坐标及接收地震波时间假定地震波在各种介质和各个方向的传播速度均相等，并且在传播过程中保持不变。现在需要根据表1中的数据确定这次地震的震中位置、震源深度以及地震发生的时间。3.1.2模型一的确定思路由于P波和S波传播的速度不同，我们采用均匀介质模型，后引入虚波，利用公式“路程=距离*时间”建立模型，提出一种线性单事件模型。3.1.3模型一的建立设震源坐标为，发震时刻为，台站i的坐标为P波、S波到时分别为，震源到台站i的距离为。采用均匀介质模型，设P波、S波速度分别为，引入虚波速度:则有:，从而（4）对站台i有：，（5）对站台j有：，（6）上面两式展开后相减，即可消去二次项，得线性方程：，（7）其中。将线性方程用于3对以上的站台，即得到一关于的线性方程组，从而反演得到震源位置。由于z:表示台站高度，当台站不具有深度时，远小于和，使得方程组的系数矩阵具有奇异性，导致z。解发散。可以通过其它方法单独求解：取震中解，将其代入(5)式求得，将各个台站求得的的平均值作为震源深度的解。这样的误差仅决定于震中的定位误差和中的到时读数误差(二者均是可控的)，从根本上解决了无法得到震源深度有效信息的问题。关于发震时刻，有方程：，（8）将方程(8)用于各个台站，求各t。的平均值作为发震时刻的解。由于完全取决于到时(见((4)式)，所以和震源的求解完全独立，其误差的唯一来源是到时读数的误差，只要到时读数足够精确，即使震源定位误差很大，也可求得很精确的。3.1.4模型一的求解采用奇异值分解(SVD)的广义逆求解。将方程组表示为:，（9）其中：，，.对矩阵A进行奇异值分解:,(10)其中，，。满足：，，奇异值是矩阵的p个非零特征值的正平方根。则方程组(9)的广义逆解为:（11）因为原始方程组即为线性，故求解过程不需要迭代，也不存在初值选取问题。并且矩阵A的计算只需要简单的算术加减，非常容易，计算量小。当A的条件数很大时，可以采用阻尼广义逆解:（12）其中I为单位矩阵，少为阻尼系数，取值随情况而变，对于震中解，取值一般在0~1之间。此时直接采用（11）式求震中解。前己述及，当台站不具有深度时，(7)式中的远小于和，这表现为矩阵A的3个奇异值中，最小值远小、于其它两个奇异值，即A的条件数很大。而对应模型空间中z方向的信息，故深度误差己将到时读数误差(包含在b中)放大了许多倍，用广义逆及阻尼广义逆均无法得到深度的有效信息。本模型充分利用台站与震源的距离可由P波和S波的到时差事先确定的特点，直接由简单的距离关系求得震源深度，避免了矩阵奇异问题。但是台站一旦具有深度，增大，矩阵A的条件数将很快减小，此时广义逆解中的深度解将获得很大改善。再次利用，可事先确定的特点，使与震源参数完全解耦，发震时刻解的精度完全不受震源定位精度的影响。3.1.5模型一的讨论选取哪两个台站的方程相减构成方程组，对定位结果有着很大的影响。将台站按照所有可能进行两两组合，这是一种操作最简单的“盲组合”，并且含有的物理信息也最丰富。实验证明这样的定位结果是最理想的图(6)。本文分别用表示震中、震源深度和发震时刻的定位误差。图6’o’表示台站,’+’表示震中。两台站间的连线表示这两个台站的方程相减所得到的方程,第4种情况含有所有可能的台站连线,即所有可能的两两组合,所得误差最小。图6台站组合图如果对一个台站的到时在其读数误差范围内进行随机扰动,得到一系列不同的到时,相当于在该台站附近增加了一系列虚拟台站,对增加后的台站进行所有可能的两两组合,(从而增加了方程个数),可以进一步提高定位结果的精度,(图7)。实际上,到时读数的确切误差并不知道,只知道一个误差范围,增加虚拟台站相当于在该误差范围内取平均效果,一般说来,这比单个的到时读数更接近准确到时,在数学形式上则表现为方程个数增加,从而所含信息量增大,故结果的精确度提高。图7误差曲线分析图考虑到表1提供的数据只有10组,我们无法利用模型一得到较精确地结果。于是我们又提出了基于几组数据即可得到较准确地数学模型。3.1.6模型二的确定思路假设采用均匀介质模型,利用空间坐标公式,求出台站到震源的距离,利用“路程=距离×时间”算出走时,再根据地震发生时间与地震台站接收时间的关系列出函数表达式。3.1.7模型二的建立假设震源三维坐标为,单位km。设发生地震的时间为2008年4月1日9时分，地震波传播速度为km/s，,i=1,2,...分别表示地震观测站点A—J的三维坐标，用表示9时分接收到地震波。根据题设条件和以上假设可得如下模型:走时=震源到观测台站的距离/地震波传播速度，观测台站接收时间=地震发生时间+走时。于是我们得到（13）其中，，，为未知参数。3.1.8模型二的求解令，，由(13)式可以看出t是x和y的函数，函数表达式中的未知参数，，，可由最小二乘法估计得到，相应的matlab程序如下：function[b,fval]=dizhenxyz=[500,3300,0300,200,0800,1600,01400,2200,01700,700,02300,2800,02900,900,03200,3100,03400,100,0];min=[21191413111410111716];sec=[9295117464714465749];time=min+sec/60;b=nlinfit(xyz,time,@dizhenmoxing,[1000100111]);functiony=dizhenmoxing(beta,x)xb=x(:,1);yb=x(:,2);zb=x(:,3);y=sqrt((xb-beta(1)).^2+(yb-beta(2)).^2+(zb-beta(3)).^2/(60*beta(4))+beta(5);运行以上程序得到：=2200.5=1399.9=35.1（14）=3=7将以上结果代入模型(13)，算得各观测站点理论上接收到地震波的时间见表2。观测站接收地震波实际时间(分)接收地震波理论时间(分)时间差(分)A21.1521.1566-0.0066B19.483319.4770.0063C14.8514.84990.0001D13.283313.27830.005E11.766711.7761-0.0049F14.783314.7881-0.0048G10.233310.23120.0021H11.766711.7677-0.001I17.9517.94640.0036J16.816716.81660.0001表2各观测站点理论上接收到地震波的时间【注:时间差=接收地震波实际时间-接收地震波理论时间】3.1.9结果的探究由上表可知各观测站接收到地震波的实际时间与理论时间之差比较小，结果(14)可以接受，即地震震中位置为原题图1中点(2200.5,1399.9)处，震源深度为35.1千米，地震发生时刻为2008年4月1日9时7分。图8为模型（13）求解后的立体曲线图。即：的曲线图图8立体曲线图附：matlab相应程序x=linspace(0,5000,50);y=linspace(0,5000,50);[xx,yy]=meshgrid(x,y);zz=7+sqrt((xx-2200.5).^2+(yy-1399.9).^2+35.1*35.1)/180;plot3(xx,yy,zz)3.2问题二的模型的建立与求解3.2.1问题的重述地震灾区的地形、地貌对抗震救灾的进展会有很大影响，根据卫星遥感和飞机航拍得到的照片可以构建灾区地形、地貌图。所构建的灾区局部地形、地貌等高线图见图9。图9灾区局部地形、地貌等高线图现在需要根据图9中的局部地形、地貌等高线图，建立数学模型，绘制出相应的三维地形、地貌曲面图。3.2.2图像的处理首先对图9进行预处理，处理效果如图10。图10灾区局部地形、地貌等高线图预处理效果图将图10保存成文件dx2.jpg，然后利用imread命令导入matlab进行自动识别，得出各条等高线所对应点的三维坐标，然后在图10中人为选定一些参照点，现在选取如下参照点：(1,1,0)(1,2,0)(1,3,0)(1,4,0)(1,5,0)(2,1,0)(3,1,0)(4,1,0)(5,1,0)(5,2,0)(5,3,0)(5,4,0)(5,5,0)(2,5,0)(3,5,0)(4,5,0)对自动识别出的各条等高线上的点和人为选定的参照点的三维坐标进行二维插值拟合，可以得到插值拟合曲面。相应的matlab程序如下：I=imread('dx2.jpg');BW=im2bw(I);[m,n]=size(BW);x_xishu=4/(n-1);y_xishu=4/(m-1);gao=[0.1:0.1:0.7];[BW,n]=bwlabel(BW);xyz=[];fori=1:n;BW1=BW;BW1(BW1~=i)=0;[y,x]=find(BW1);xyz=[xyz;110;120;130;140;150;210;310;410;510;520;530;540;550;250;350;450];x=xyz(:1);y=xyz(:2);z=xyz(:3);[x1,y1]=meshgrid(1:0.1:5);z1=griddata(x,y,z,x1,y1,'v4');figuresurf(x1,y1,z1)zlim([00.8])插值拟合曲面，即与图9中的局部地形、地貌等高线图相应的三维地形、地貌曲面图如下：图10三维地形、地貌曲面图3.3问题三的模型的建立与求解3.3.1问题的重述现在需要查阅有关资料，了解地震波在各种介质和各个方向的传播速度问题，给出合理假设，根据表1中的数据确定这次地震的震中位置、震源深度以及地震发生的时间。3.3.2模型的确定思路在假设模型的过程中，我们查阅了大量的关于地震问题的文献，我们发现Geiger模型对我们的问题具有指导性意义。我们将问题一的解作为Geiger方法的初值，因为初值离准确解已经相当近，所以能够保证迭代过程收敛至全局极小值，有效地提高了定位精度。3.3.3模型的建立设个台站的观测到时为，求震源及发震时刻使得目标函数（15）最小。其中为到时残差（16）为震源到第i个台站的计算走时。3.3.4模型的求解使目标函数取极小值也即（17）其中。为方便，记（18）则由(18)式,在真解附近任意试探解及其校正矢量满足(19)也即（20）由的定义可得公式(20)的具体表达式（21）若偏离真解不大，则和较小，可忽略二阶导数项，（21）式被简化为线性最小二乘解：（22）以矩阵形式表示，上式为（23）其中，.若二阶导数项不可忽略，则(12)式给出非线性最小二乘解（24）通常各台站的到时数据具有不同的精度，如果不加以区别，则具有较低精度的数据将严重干扰结果的精度，这一问题可以通过引入加权目标函数来解决。设各台站到时残差的方差为，引入加权目标函数，（25）按照上述同样的步骤，通过求(25)式的极小值，得到如下加权线性最小二乘解，（26）其中为加权方差矩阵：。由方程(24),(25),或(26)求得后，以作为新的尝试点，再求解相应方程。如此反复迭代，直至或足够小（或满足一定的循环结束条件），此时即得估计解。然而，Geiger方法也有它的不足之处。Geiger方法强烈地依赖于初值的选取，如果初值不够好，会导致迭代过程的发散，根本无法求得正确解，当初值离准确解较近时，能保证迭代收敛，且迭代次数较少。所以本题我们选取了模型二的解值作为我们的初值。3.3.5模型的进一步探究我们想到能否用另一个模型可以补充Geiger方法的不足之处，考虑到地震波在各方向传播的问题，我们采用震源位置与速度结构的联合反演法（SSH），SSH法不需要对波速进行校准，同时还可以获得有关速度结构的很多信息，是目前被广泛使用的一