地震波对房屋晃动的动态模拟

摘要：地震波对房屋晃动的动态模拟主要的目的是用于教科书上的p波和s波的科普，当p波和s波到达地面时，地面会分别出现不同的晃动而做出的研究，为了形象的模拟地面的晃动情况，就用模拟房屋的晃动来表示地面的情况。首先通过读取地震台站的数据，然后根据震源，震中距以及震源深度，来计算地震波的p波和s波的传播速度以及走时。最后模拟出房屋，使其房屋根据p波和s波的到达先后而出现不同的晃动。

关键词：地震波；震源；震中距；震源深度；传播速度；走时

中图分类号：TP311 文献标识码：A 文章编号：1009-3044（2018）26-0219-02

在《地震概论》这本书中我们有学到p波和s波，书中是按画图的形似描绘出p波和s波的性质特征，虽然描述的很详细，但是静态的图片往往不如动态画面更加清楚、更加直观的体现p波和s波的震动方向和传播速度。因为现在的时代是一个快节奏时代，图片和文字以及静态的描述一个事物往往更不容易使人们去接受它。所以地震波对房屋晃动的动态模拟可以让人们更加容易的接受和理解p波和s波。

通过地震波对房屋晃动的动态模拟，对于授课老师说，可以以这次模拟作为辅助材料，以演示程序的方式给同学们来讲解p波和s波，使老师讲解起来更加轻松，同学们也能更好的理解p波和s波。在宣传地震波知识和传播地震波知识也都可以使用这次模拟来给人们讲解。

地震波按传播方式分为三种类型：纵波、横波和面波。纵波是推进波，地壳中传播的速度大约为5.5～7千米/秒，最先达到震中，又称p波。横波是剪切波，在地壳中的传播速度大约为3.2～4.0千米/秒，第二个达到震中，又称为S波。面波又称为L波，是由纵波与横波在地表相遇后激发产生的混合波。

地震波是一种机械运动的传播，产生于地球介质的弹性。它的性质和声波很接近，但普通的声波在流体中传播，而地震波是在地球介质中传播，所以要复杂的多。但这次模拟主要是研究S波和P波的性质，并且，地球介质具有多元性和复杂性，所以抽取最主要的部分进行简化，即不考虑P波和S波在传播过程中由于地球介质而出现的弹性波。在计算地震波的时候和光波有些相似之处。波动光学在短波的情况下可以多渡到几何光学，从而简化了计算；同样地，在一定条件下地震波的概念可以用地震射线来代替而形成了几何地震学。

这次模拟首先是获取地震波的数据，对其进行数据的转换，分别计算出P波和S波的传播速度。然后对其地壳进行分层，再根据震中距计算入射角等参数。根据射线在地下的传播路径来描述p波和s波，当射线达到房屋时，通过与房屋的衔接，可以看到房屋晃动的模拟。

1 数据说明

程序的数据从国家地震科学数据共享中心离线获取，地震波形是2015年11月28日2：10：37.78时刻发生在河北省唐山市丰南区Ms3.4级地震为例，其震源深度为21km，选取的台站为LUT台、BET台。

2 地震波数据的转换

离线下载的数据格式不能用于MATLAB软件，所以我们用MSDP软件打开地震数据，将其导出成Ascii码，文件转成为rar格式。然后再解压，即得到每个台站各个通道的txt格式数据，在导出波形数据时，从地震发生时刻开始导出。由于波形数据分为三个方向，南北，东西，垂直向，程序中我们选用的通道是垂直向的，也就是说在动画中出现的波形都为垂直向通道的，且都是以直达波形形式表示射线在平行层中传播。

3 地震波速度的确定和对地壳进行分层

v=[2.5，5.3，6.1，6.6，7.2]； %地壳模型的各层速度

deep=[0，1，3，24，38，46]； %地壳模型的各层厚度

for ii=1：length（v）

%两个界面所夹的层为均匀层具有相同的速度

deep1=[deep1，deep（ii），deep（ii+1）]；

%在绘图的时候根据颜色来区分不同的层级 vcolor=[vcolor；ones（1，nlength）*v（ii）；ones（1，nlength）*v（ii）]；

end

h0=38.1； % 地震深度（为了满足需要，我们对震源深度进行了修改）

Dz = h0；

SLayerNum = 1；

for l=2：length（deep）

%如震源深度小于速度模型中某层，那么震源在该层的上一层

if （h0 < deep（l））

Dz = h0 - deep（l-1）；

%计算震源所在层

SLayerNum = l-1；

break；

end

end %如果震源深度大于速度模型中最深层，计算震源距层顶距离和所在层

if h0 > deep（length（deep））

Dz = h0 - deep（length（deep））；

SLayerNum = length（deep）；

End

4 地震波发射的角度

x=30； %设置震中距为30

n=n+1；

%计算迭代离源角的最大值

Maxthetaj = atan（x/Dz）；

%計算迭代离源角的最小角

Minthetaj = atan（x/h0）；

thetaj= （Maxthetaj + Minthetaj）/2； %计算入射角初值

%根据新的入射角和速度模型，重新计算震中距

EpiDis=Dz*tan（thetaj）；

for l=SLayerNum-1：-1：1

tanThetajl = sin（thetaj）/sqrt（（v（SLayerNum）/v（l））^2 - （sin（thetaj））^2）；

EpiDis= EpiDis + （deep（l+1）-deep（l））*tanThetajl；

end

%迭代求出入射角（逼近到系统设置值迭代结束）

while （abs（EpiDis - x） > 1.0e-10）

%根据震中距计算入射角

if x > EpiDis

Minthetaj = thetaj；

elseif x < EpiDis

Maxthetaj = thetaj；

end

thetaj = （Maxthetaj + Minthetaj）/2；

%根据新的入射角，重新计算震中距

EpiDis=Dz*tan（thetaj）；

for l=SLayerNum-1：-1：1

tanThetajl = sin（thetaj）/sqrt（（v（SLayerNum）/v（l））^2 - （sin（thetaj））^2）；

EpiDis= EpiDis + （deep（l+1）-deep（l））*tanThetajl；

end

%如大角与小角差值小于10E-10，迭代结束

if（abs（Minthetaj - Maxthetaj）<10e-10）

break； %退出迭代

end

end %WHILE循环结束

经过震源，震中距和最右离源角等参数的确定，再进行了一系列的图像处理以后，就产生了以下的效果：

这是一个二维坐标图，其中最底下是指用颜色来表示地殼的厚度，每一层都是用了不同的颜色来区分，这样能使地壳分层更加的直观。右下角是图例，蓝色代表的是p波路径，绿色的代表的是s波路径，在前面我们已经给p波和s波的不同层级的地壳速度的计算，程序运行的结果为：先是p波到达地面，此时房屋开始上下晃动。然后紧接着是s波的达到，由于p波的作用还在，所以p波和s波开始杂糅在一起，使房屋开始了无规律的剧烈晃动，最后，p波的作用慢慢褪去，由于房屋还有s波的作用所以房屋会出现左右晃动。

由程序结果可以知道：p波的速度比s波的速度，所以p波会先达到地面。p波的震动方向是上下方向，s波的震动方向是左右的。往往上下晃动的时候对房屋是没有很大的影响的。所以在s波达到之前，是最佳逃生的时间。

由于数据等一些原因，房屋和地壳都是以二维的效果简单的表示了一下，在视觉上没有像三维的效果图一样的那么生动。并且没有考虑到房屋的抗震系数也就没设计出房屋倒塌的情景，所以只适合在宣传地震波知识和在地震发生时我们该怎么应对的时候作为辅助的课件演示。

在自然灾害当中，其中地震的杀伤力最大，因为地震的不确定性和传播快的特点。所以希望这次的模拟能给每一个人都能够区分p波和s波，并且遇到了地震的话，能够在最佳的逃生时间内自救。

参考文献：

[1] 赵晓燕，于仁宝.地震概论[M].北京：清华大学出版社，2013：1-177.

[2] 万永革.地震学导论[M].北京：科学出版社，2016：1-469.

[通联编辑：唐一东]