地震震中位置的确定和三维地形图的绘制

地震震源位置确定及三维地形图的探讨
摘要
本文根据多个地震观测点接收地震波的时间，在假设条件下，对地震的震中位置、震源深度及地震发生的时间和三维地形图的绘制进行了建模求解。

并在问题（3）中，增加了一些实际考虑因素，将模型进行拓展，得到一个更合乎实际的模型。

针对问题一，我们以G 点为参考点，对题目中所给的数据进行预处理，运用平面解析几何的数学方法，推导出各个观测点到震中位置的关系式。

据此，可以建立确定震中位置、震源深度及地震发生时间的数学模型：
()()222)(h y y x x t t v i i i +-+-=
+⋅
针对问题二，利用MATLAB 软件中的imread 函数将等高线图中的全部数据（以像素为单位）存储到三维色彩分量矩阵中，借助MATLAB 的矩阵处理功能，对等高线颜色的数据分离处理，得到色彩分析表（见表3）。

随后，寻找等高线与其对应地理坐标的关系，建立起等高线图和三维地形图坐标变换的数学模型：
],[i i y x i E →
→
==]183
2
,771[i y x i →
→)7,3,2,1( =i 从而求得高度不同的等高线对应的地理坐标。

最后，借助MATLAB 软件的
griddata 函数cublic 插值，拟合出局部地形等高线图的三维地形图（见图2）。

针对问题三，我们查阅了相关资料，发现纵波传播速度最快，最容易被地震观测点观测到。

因此，本文只考虑纵波对震中等位置确定的影响，提出了 5.3.1中的五个假设。

并根据假设，对模型一进行改进，可建立地震的震中位置、震源深度及地震发生的时间确定的数学模型三：
()221222
2)()()3
2(5.522
h h y y x x h t t h i i i ++-+-=
-+⋅+
运用MATLAB 软件可以解得：震中位置坐标为（2244.2759，1134.3445）千米；地下水层的厚度404.243千米，花岗岩的厚度943.235千米，震源的深度1347.478千米；地震发生的时间为：2011年4月1日9时7分20秒。

关键词：地震震中位置的确定 三维地形图 坐标变换 cublic 插值
一、问题重述
1.1背景
2008年5月12日发生在四川汶川的一场大地震给我们每一位中国人带来了巨大的伤痛，痛定思痛，我们应该为减少震后灾害做些事情。

当地震发生时，震中位置的快速确定对第一时间展开抗震救灾起到非常重要的作用，而震中位置可以通过多个地震观测站点接收到地震波的时间推算得到。

这里假定地面是一个平面，在这个平面上建立坐标系见图1。

图1中给出了10个地震观测站点（A —J ）的坐标位置。

500
1000
15002000
2500
3000
3500
0500
1000
1500
2000
2500
3000
3500
X （单位：千米）
Y （单位：千米
）
图1 地震观测站点示意图
2008年4月1日某时在某一地点发生了一次地震，图1中10个地震观测站点均接收到了地震波，观测数据见表1。

地震灾区的地形、地貌对抗震救灾的进展会有很大影响，根据卫星遥感和飞机航拍得到的照片可以构建灾区地形、地貌图。

所构建的灾区局部地形、地貌等高线图见图2。

注：（1）图2中等高线标注数字的单位为千米；
（2）图1和图2不属于同一坐标系。

X （单位：千米）
Y （单位：千米）
1 1.52
2.53
3.54
4.55
图2 灾区局部地形、地貌等高线图
1.2地震观测站坐标及接收地震波时间的数据
表1 地震观测站坐标及接收地震波时间
1.3要解决的问题
问题一：假定地震波在各种介质和各个方向的传播速度均相等，并且在传播过程中保持不变。

请你根据表1中的数据确定这次地震的震中位置、震源深度以及地震发生的时间（不考虑时区因素，建议时间以分为单位）。

问题二：请你根据图2中的局部地形、地貌等高线图，建立数学模型，绘制出相应的三维地形、地貌曲面图。

问题三：查阅有关资料，了解地震波在各种介质和各个方向的传播速度问题，给出合理假设，根据表1中的数据确定这次地震的震中位置、震源深度以及地震发生的时间
二问题分析
2.1对问题（1）的分析
问题一要求我们在假设地震波传播速度不变时，根据附表中的数据，确定这次地震的震中位置、震源深度及地震发生的时间。

首先，根据题中所给的地震观测图，假设震中的坐标为O(x，y)，震源深度为h千米，G点接收到地震波的时间为地震发生后t秒，且地震波的传播速度为v千米/秒。

其次，以G点为参考点，对题目中所给的数据进行预处理，得到表1。

运用平面解析几何的数学方法，做出地震波传播到某观测点的简化几何图，推导出各个观测点到震中位置的关系式，可以建立震中位置、震源深度及地震发生时间的数学模型。

最后，任意选取五个地震观测点（E,H,C,D,G），将其坐标代入模型一中，运用MATLAB软件，求出震中位置、震源深度及地震发生的时间。

2.2对问题（2）的分析
问题二要求我们根据图2中的局部地形、地貌等高线图，建立数学模型，绘制出相应的三维地形、地貌曲面图。

首先，我们利用MATLAB软件中的imread函数将等高线图中的全部数据（以像素为单位）存储到三维色彩分量矩阵中，借助MATLAB的矩阵处理功能，按照等高线颜色的不同，进行数据的分离处理。

其次，通过寻找等高线与其对应地理坐标的关系，建立起坐标变换的数学模型，从而求得高度不同的等高线对应的地理坐标。

最后，提取并分离得到的数据点后，借助MATLAB软件的griddata函数进行cublic插值，拟合出局部地形地貌等高线图的三维地形图。

2.3对问题（3）的分析
问题三要求我们根据查阅的相关资料，给出合理的假设，并根据表1中的数据确定这次地震的震中位置、震源深度以及地震发生的时间。

首先，我们从百度百科、《现代地震学教程》和《地震勘测原理》的书中了解到：纵波在地下水层中的传播速度为1.5千米/秒,在花岗岩中的传播速度为5.5千米/秒，且纵波在地壳中的传播速度远大于横波。

由于纵波传播速度最快，最容易被地震观测点观测到，因此，只考虑纵波对震中等位置确定的影响。

其次，根据查的资料和4.2的名词解释，我们做出各个地震观测点记录的观测时间均为接受到纵波的时间，及从地震波传播的方向，传播时经过的介质（为均匀介质）以及传播的速度等做了具体合理的假设。

最后，画出简化分析图，运用解析几何的方法，可以推导模型三，并运用MATLAB软件对模型三进行求解，得到结论。

三、模型假设
（1）假设所有观测点在同一平面上；
（2）假设地震波在各种介质和各个方向上的传播速度保持不变；
（3）假设各个观测点所测量的数据准确无误且地形在地震中没有剧烈的变化；
（4）相邻观测点之间的高程变化是线性的，可以利用差值求取两个点之间的点；
（5）三次样条插值时，各个区间段的函数是连续的。

四、模型准备
4.2名词解释：
1、地震波的传播方式：
地震波按传播的方式分为三种类型：纵波、横波和面波。

纵波为推进波，地壳中传播速度为5.5~7千米/秒。

最先到达震中，又称为P波，它使得地面发生上下振动破坏性叫弱。

横波是剪切波：在地壳中的传播速度为
3.2~
4.0千米、秒，第二个到达震中，又称为S 波，它使得地面发生前后、
左右抖动，破坏性较强。

面波又称为L波，是由纵波与横波在地表相遇后激发产生的混合波。

其波长大、振幅强，只能沿着地表面传播，是造成建筑物
剧烈破坏的主要因素。

2、等高线：
等高线指的是地形图上高程相等的各点所连成的闭合曲线。

把地面上海
拔高度相同的点连成的闭合曲线。

垂直投影到一个标准面上，并按比例缩小画在图纸上，就得到等高线。

等高线可以看作是不同海拔高度的水平面与实际地面的交线，所以等高线是闭合的曲线。

在等高线上标注的数字为等高线的海拔高度。

五、模型建立及求解
5.1震中位置、震源深度及地震发生时间的确定
5.1.1震中位置、震源深度及地震发生时间数学模型的建立（模型一）
根据2.1的问题分析，我们可以知道，假定地震波在各种介质和各个方向的传播速度相等，且在传播过程中保持不变。

据此，设震中的坐标为O ),(y x ,震源深度为h 千米，G 点收到地震波的时间为地震发生后的t 秒，地震波传播的速度为v 千米/秒，地震波从震源（O ）出发传播到某个地震观测点S ),(i i y x ''所需要的时间比到G 点的时间晚i t 秒。

其次，对题目中的数据进行预处理：以G 点为参考点，将题中所给表的时间换为以秒为单位，可以得到表2：
由此，可以做出地震波从波源传播到某观测点的简化几何图如下：
图1 地震波从波源传播到某观测点的简化几何图
根据图中的关系，地震波从O 点传播到S 点的总路程为OS,运用解析几何的方法，可以推导出下面的公式（模型一）：
()()222)(h y y x x t t v i i i +'-+'-=
+⋅
5.1.2震中位置、震源深度及地震发生时间数学模型的求解
借助题目中所给的数据，选取E,H,G,D,C 五个观测点作为参考点（任取），将这五个观测点的坐标代入模型一，可得：
对E 点：()()2
22)700(170092h y x t v +-+-=+⋅
对H 点：222)900()2900()92(h y x t v +-+-=+⋅
对G 点：222)1900()2500(h y x t v +-+-=
⋅
对D 点: 222)2200()1400()183(h y x t v +-+-=+⋅ 对C 点：222)1600()800()277(h y x t v +-+-=
+⋅
在MATLAB 软件中运行附录中的程序一，可以解得：震中位置、震源深度及地震发生的时间，如表3：
5.2三维地形、地貌曲面图的绘制
5.2.1等高线图的各处等高线处数据的导入与提取
先将等高线图转换成bmp 格式，保存在MATLAB 的当前目录下，利用MATLAB 软件中的imread 函数将等高线图中的全部数据（以像素为单位）存储到三维色彩分量矩阵。

（由于bmp 格式的图片以像素的个数和色彩的深度确定的，本文以宽度422像素，高度353像素，色彩深度为24位的图片为例。

）
因此，导入的三维色彩分量矩阵为：353*422*3，共存储了353*422个像素
点的色彩数据。

利用matlab 图形窗口中的Data Cursor 工具可以读取图像上各点的数据，即得到了各等高线处的色彩分量（见表4）。

5.2.2等高线图与三维地形图间坐标变换的数学模型（模型二）
先利用MATLAB 软件内置的find 函数将导入的矩阵拆成三个二维的矩阵I1，I2，I3分别代表三个色彩的分量。

1、为了更好地反映图片中各高度的等高线上样本点的矩阵位置，我们用find 函数找到对应的矩阵位置后，将其返回到三个列数为2的矩阵中，取其交集，得到列数为2的矩阵：
],[i
i y x i E →
→
=)7,3,2,1( =i 2、为了使矩阵的位置变为实际意义的地理坐标，我们利用原图中单位像素对应的实际距离进行坐标变换，利用matlab 图形窗口中的Data Cursor 工具，提取图中坐标轴顶点的位置（见表5）
320-12=308，所以单位水平像素所对应的实际长度为77
1
km ；垂直坐标4km 对应
的像素差为408-42=366，所以单位垂直像素所对应的实际长度为183
2
km 。

3、等高线图与三维地形图间坐标变换的数学模型
将2中的变换系数与式(1)相结合，就可得到各等高线上的像素点的地理坐标 ：
],[i i y x i E →→==]183
2
,771[i y x i →
→)7,3,2,1( =i 由上面的模型，可以得到绘制三维地形图所需要的两个向量（样本点地理位
置的横纵坐标）为：
127127[,,,]
[,,,]T T T T T T T T x x x x y y y y ⎧=⎪⎨
⎪=⎩
5.2.3三维地形、地貌曲面图的绘制
根据5.2.1和5.2.2的数据处理和模型建立，结合调用MATLAB 软件中的散点插值函数griddata ，我们编写出附录中三维地形地貌曲面图的绘制程序。

将其在matlab 软件中运行后，再MATLAB 的命令窗口输入：[X1,Y1,Z1]=griddata(X,Y,Z,(linspace(1,5))',linspace(1,5),'cubic'); mesh(X1,Y1,Z1)，可以得出该局部地区的三维地形、地貌曲面图（见图2）
图2 该局部地区的三维地形、地貌曲面图
5.3结合假设条件，对震中位置、震源深度及地震发生时间进行确定
5.3.1震中位置、震源深度及地震发生时间数学模型的求解
我们从百度百科、《现代地震学教程》和《地震勘测原理》的书中了解到：纵波在地下水层中的传播速度为1.5千米/秒,在花岗岩中的传播速度为5.5千米/秒，且纵波在地壳中的传播速度远大于横波。

因此，在建立模型三时，先假设地震观测点记录的观测时间为接受到纵波的时间。

结合4.2的名词解释，可以做出如下假设：
（1）假设纵波地震波从震源传播到某个观测点都需要经过两层介质（地下水层和花岗岩层，且先经过水层比花岗岩层更接近震源）且两层介质均是均匀的，不影响地震波的传播速度；
（2）假设从震源到地平面，地下水层与花岗岩层的厚度比为：7:3:21 h h ； （3）假设地震波穿过地下水层和花岗岩层的时间分别为：21,t t ；
（4）假设地震波穿越两层介质到达地面的地震观测点是沿直线传播的，在两层之间发生的速度变化忽略不计；
（5）假设地震波传播方向为：其到所有观测点的直线与地平面所成的角度都为o 45=β。

根据上面的假设，可以做出假设的简化分析图，以震源O 到达某个地震观测点S1为例：（见图3）
图3 地震波从震源历经两层介质传播到某观测点的简化几何图
根据图中的关系，地震波从O 点传播到S1点的总路程为OS1,运用解析几何的方法，可以推导出下面的公式：
1122
11
21212cos 45t :3:7
i v t v t v t h t t h h ︒
⎧⋅+⋅=⎪⎪⋅=⎨⎪+=⎪
=⎩ 由上面的公式，我们可以推导出模型三为：
()221222
2)()()3
2(5.52
2
h h y y x x h t t h i i i ++-+-=
-+⋅+
5.3.2震中位置、震源深度及地震发生时间数学模型的求解
借助题目中所给的数据，选取I,H,G,D,C 五个观测点作为参考点（任取），将这五个观测点的坐标代入模型三，可得：
对I 点:()2212222)()3100(3200)3
292(5.52
2
h h y x h t h ++-+-=-+⋅+
对H 点：221222
2)()900()2900()3
2463(5.52
2
h h y x h t h ++-+-=-+⋅+
对G 点：2212222
)()1900()2500()32(5.522h h y x h t h ++-+-=-⋅+ 对D 点: 2212222
)()2200()1400()32183(5.522h h y x h t h ++-+-=-+⋅+ 对C 点：221222
2)()1600()800()3
2277(5.522h h y x h t h ++-+-=-+⋅+ 在MATLAB 软件中运行附录中的程序三，可以解得：震中位置、震源深度及地震发生的时间，如表6：
并绘制出震中的位置观测图（见图3）
050010001500
20002500300035000500
1000
1500
2000
2500
3000
3500
x y
图3 震中的位置观测图
图中，红点处为震中的位置，其坐标为(2244.2759,1134.3445)千米。

六 模型评价
本文建立了三个模型对地震的震中位置、震源深度及地震发生的时间进行了逐步深入的分析求解，并绘制了局部地区的三维地形地貌图。

对于模型一，我们考虑通过简单的几何关系，建立起确定震中位置等问题的模型。

该模型的优点是
简单易懂，对数据的处理直观明了，较为准确地得到所要求得数据。

其不足在于，没有考虑其它外界因素对地震波传播的影响，有一定的局限性。

对于模型二，所建立的等高线图和三维地形图的坐标变换模型，能够巧妙地从像素的数据出发，抽取合适的样本点最终得到较为满意的三维地形图。

该模型具有广泛的应用性。

对于模型三，我们只考虑纵波对震中等位置确定的影响，对问题缩小来处理，使得模型得到的结果较为精确。

但是，考虑实际的因素过少，实际意义不大。

因此，在模型三的改进中，我们提出了在模型三的基础上，增加横波和介质分层增多等因素，使得模型的更具有实际意义，预测的位置更为精确。

七、参考文献
[1]邢永强、李金荣、张天义、杨振放、潘涛，地震预报方法研究[J]，安徽农业科学，12845—12846，2008年；
[2]刘静，利用MATLAB由山区等高线图模拟三维地形图，中国科技论文在线，,2011年8月31日；
[3]陆基孟，地震勘探原理，山东：中国石油大学出版社，2009年3月；
[4]刘卫国，MATLAB程序设计教程（第二版），北京：中国水利水电出版社，2010年2月；
[5]胡蓉，MATLAB软件与教学实验，北京：经济科学出版社，2010年8月；
[6]百度百科，地震波及等高线名词定义的解释，，2011年8月31日；
[7]周仕勇、许忠淮，现代地震学教程，北京：北京大学出版社，2010年12月1日。

八、附录。