用地震台站记录得到区域地壳厚度与横波速度

这几天在上大三暑假的小学期，一共上了三位老师的课程，分别是大地电磁、地震勘探和天然地震的，各有千秋（成语用的不是很好），让我大开眼界。趁着记忆新鲜，先把天然地震涉及到的课题记录下来，以备后日查看，也能成为不止保存在我回忆中的印刻。大地电磁有一部分是用python写一个读取某类型数据文件的小程序，我记得我写的非常开心，但是最后许多功能其实还没成功实现，所以之后可能会再抽一些时间去研究。
这篇文章主要是阐述一下如何通过同一个台站接收记录到的三分量地震数据得到比较可靠的地壳厚度与横波速度（纵波速度假定6.1km/s，计算过程中其实代进去的是波速比k=Vp/Vs）。涉及到的地震软件有Seismic Analysis Code(SAC)、TauP以及Computer Programs in Seismology(CPS)，编程和处理都是在我Mac上进行的，Linux上也完全no普拉卜冷！
然后要说啥呢，我想想。嗷，先概述一下原理，再列个目录。
基于以下三个计算震相走时的公式：

H-地壳厚度 k-波速比 p-射线参数

通过遍历可能的H和k值，分别计算每一种情况下的三个到时，用这三个到时去对应每个Event处理出来的接受函数得到振幅。倘若H和k值接近真实值，那么振幅之和应当是极大值甚至是最大值。最后绘制H、k、振幅之和的等值线图，就可以直观的观察到不同H和k的效果。因此，大体上按照以下几个步骤进行：

  • 对原始三分量地震记录进行处理，得到接收函数
  • 凭借震中距与震源深度得到射线参数
  • python计算绘等值线图

对原始三分量地震记录进行处理得到接受函数

说到这个三分量地震记录，首先我们得知道数据文件的格式，一般由道头数据和地震数据组成，跟segy的款式有点像，但是应该是不完全一样的。具体的描述在Seisman的SAC参考手册上都有，在这里就不班门弄斧。总之，每个地震Event都拥有一个独立的目录，该目录下就是三个位移分量数据。在XQuartz上进入某Event目录后，直接‘sac’调用SAC软件（会出现一段软件的版本和开发者信息），就可以读入你要的数据文件，并且绘制出时间序列。

bash-3.2$ sac
 SEISMIC ANALYSIS CODE [11/11/2013 (Version 101.6a)]
 Copyright 1995 Regents of the University of California

SAC> r *.SAC
BHE.SAC BHN.SAC BHZ.SAC
SAC> qdp off
SAC> p1

此处命令较为常用，都采用了缩写，r就是read--读入文件，p1就是plot1--把所有的时间序列一并画出来，qdp off是关闭快速绘图的功能，具体的我没有多看，但从off前后的图像对比我猜测qdp的功能就是为了快速出图所以略去一些高频成分。得到了这样的图像：


原始三分量.png 横坐标时间 纵坐标振幅

这里的数据是我已经rmean过的，因为地震仪质心会不断偏移，所以时间序列的中心轴不在0振幅上是非常常见的，为了校正可以在SAC里调用rmean命令去除均值。另外在第三张图中有两条标记线P和S，这是在终端中调用TauP函数得到的震相初至标记，等会儿会说到。
现在需要对三分量中的水平数据做旋转，把正南北、正东西的数据转换为径向与切向的数据，即做坐标系的变换，使用的是SAC中的rotate命令。

SAC> r BHE.SAC BHN.SAC
BHE.SAC BHN.SAC
SAC> rotate to gcp
径向与切向时间序列.png

可以明显看到径向的直达P波振幅较大，两个方向的直达S波都比较强，径向上的甚至更强一些。我记得我去问老师：为什么径向上也会有这么强的直达S波。后来我知道径向上的主要是SV波的分量，切向上主要是SH波。
我们的目的是得到径向方向上的接收函数，所以要用垂直向对径向做反褶积。但在此之前我们先把垂向和径向序列中直达P波段的时间序列截取出来，用短一点的时间序列来做反褶积。这里就要用到P波到时作为参考时间。那么怎么才能知道P波是什么时候到的呢？当然可以肉眼提取（准确度和瞎眼度另谈），不过我们这里采用TauP中的命令来自动标记。

$ taup_setsac -evdpkm -ph P-1,S-2 BHR.SAC

BHR.SAC是径向文件（Radial），BHT.SAC是切向文件（Tangential）（这里没用）
-ph P-1,S-2是选择你想要哪个震相的到时标记，这里让P波到时标记为T1MARKER，S波到时标记为T2MARKER，放入道头信息中（其实只需要P波到时就可以，但是多写一个有对比好理解）。
有了P波到时，就能够截取出P波前20s后100s段的时间序列：

SAC> cut t1 -20 100
SAC> w wR wZ

w即write，将切完的内存中的数据写入wR和wZ文件。
接下来！终于可以进行反褶积了！（呼，原来写这种教程这么慢，给我写饿了）反褶积用啥呢？用CPS的函数：

$ saciterd -FN wR -FD wZ -N 400 -D 10 -E 0.0001 -ALP 2

-FN指作为反褶积分子的数据文件，-FD是分母，-N是最多迭代次数，-D是时间序列预留空出来的时间（我不确定是两边还是光左端），-E是迭代要求误差大小（我后来就没一个满足这个要求的），-ALP是与高斯滤波器有关的选项。得到了径向的接收函数（生成的一系列文件中那个decon.out）：
接收函数.png （0时刻是PP波的到时）

因为有三十几个Event，不想一个一个读到SAC里面一个一个处理，就写了个脚本：

#!/bin/bash
SAC_DISPLAY_COPYRIGHT=0

echo "EVENTS:" > ~/Desktop/lists.txt
path=~/Desktop/HAZ.data3
cd $path
for eve in `ls`; do
    cd $path/$eve
    echo "$eve" >> ~/Desktop/lists.txt
    sac << EOF
      read *.SAC
      rmean
      write over
      read *BHE* *BHN*
      rotate to gcp
      write Radial.SAC Tangential.SAC
      quit
EOF
    taup_setsac -evdpkm -ph P-1,S-2 Radial.SAC Tangential.SAC *BHZ*
    sac << EOF
      read Radial.SAC Tangential.SAC *BHZ*
      cut t1 -20 100
      write wR wT wZ
      read wR wT wZ
      bp c 0.02 5 n 2 p 2
      write over
      quit
EOF
    cd $path/$eve
    mkdir wR_RF
    cd $path/$eve/wR_RF
    saciterd -FN $path/$eve/wR -FD $path/$eve/wZ -N 400 -D 10 -E 0.0001 -ALP 2
done

这里有个坑，sac << EOF下一行要再空两格，不然运行就会说找不到这些命令，我找这个错误找了好久，找到之后又给我气了好久。
最后在三十几个Event里选了9个反褶积效果较好的参与后续活动，其他Event一人一个谢谢参与奖。

凭借震中距与震源深度得到射线参数

这一部分其实不复杂，就是通过TauP输入每个Event的震中距和深度，然后把TauP计算好的射线参数拿出来，以供之后算t备用。手动输入输出输入输出，最多五分钟就能完事儿。不过有位老师跟我说过他的一个缺点：做啥数据处理都想用批处理，经常把5分钟的工作量硬生生拖成一个小时甚至半天。结果把我也传染了，多于5个的重复处理就非得用shell做个流程（微笑）。

#!/bin/bash
SAC_DISPLAY_COPYRIGHT=0 #关闭调用sac就会出现的提示信息，免得眼花
path=~/Desktop/data  #data目录里是挑选出来的Event，个个皆是人中龙凤
cd $path
for eve in `ls`; do
    cd $path/$eve
    sac << EOF
      read ./Radial.SAC
      taup_time -ph P -h &./Radial.SAC,EVDP& -km &./Radial.SAC,DIST& --rayp -o p_out
      quit
EOF
    #echo "p_list:" > ~/Desktop/p_list.txt
    #echo `pwd`
    echo "$eve  `awk '{printf $1}' p_out.gmt`" >> ~/Desktop/rayp_list.txt
done

这个taup_time就是用来计算射线参数等数据的函数，-ph P就是只做P波震相的射线参数，-h是震源深度，-km是震中距，--rayp是只输出射线参数，-o是输出到哪个文件里。
这里值得提的一点是，我一开始折腾了好久怎么把SAC的道头变量拿到外面来，这样就能传入taup_time做计算，但是当时就是死活找不到方法（倒也没有寻死觅活的）。后来福至心灵想既然拿不到外面来能不能在SAC里面做taup_time，发现是可以的（但是同学用Linux说不行，不知个中真相），所以就直接用&file,varname&拿到了深度和震中距值。但是！！！就在几个小时之后，在Seisman的参考手册上看到了我想要的，用saclst工具就可。但我的心已经被伤了，再挽留也回不来了。这玩意儿以后要用以后再说吧，👴找着比你更听话的了。

python计算绘等值线图

第三站 大制作（估价10元）
直接放代码（在jupyter notebook中编译）

filepath = '/Users/husir/Desktop/rayp_list.txt'
rpl_txt = open(filepath,'r')
rpl = rpl_txt.readlines()
for p in rpl:
    print(p)

for i in range(len(rpl)):
    rpl[i] = rpl[i].split()
    rpl[i][1] = float(rpl[i][1]) # 单位是s/degree
    rpl[i][1] = rpl[i][1]/111.1849 # 单位是s/km
for p in rpl:
    print(p)

import os
import numpy as np
path = "/Users/husir/Desktop/data"
events = os.listdir(path)
del events[1] # 删除.DS_Store这个目录
#print(events)
Vp = 6.1 # km/s
H = np.arange(24,56,2) # 地壳厚度
k = np.arange(1.4,2.4,0.1) # Vp/Vs波速比
Amp_sum = [] # 最后把指定H&k下所有Event的所有震相到时t在接受函数中对应的振幅都加起来\
             #储存在这里，作为衡量H和k是否准确的标准，即后面等值线图的Z（这是一个二维数组）

for iH in H:
    Amp_sum.append([])
    for ik in k:
        A_sum = 0
        for ev in events:
            for i in range(len(rpl)):
                if rpl[i][0] == ev: # 防止导入的Event顺序跟rayp不同
                    rp = rpl[i][1]
                    break
            # 计算三个震相的到时（相对于PP波的到时）：tPS, tPpPs, tPsPs+PpSs
            tPS = iH * (np.sqrt(ik**2/Vp**2-rp**2) - np.sqrt(1/Vp**2-rp**2))
            tPpPs = iH * (np.sqrt(ik**2/Vp**2-rp**2) + np.sqrt(1/Vp**2-rp**2))
            tPsPs_PpSs = 2 * iH * np.sqrt(ik**2/Vp**2-rp**2)
            #print(tPS)
            
            RF_path = os.path.join(path,ev,"wR_RF/rf2.0") # 到这个接受函数文件中提取相应震相计算到时对应的振幅
            rfp = open(RF_path,'r') # 从道头中可以看到，起始时间是-10s，结束时间是153.83s，采样率是0.01s，固共有16384个数据
            rfl = rfp.read()
            rfl = list(map(float,rfl.split()))
            # 得到APS, APpPs, APsPs+PpSs
            APS = abs(rfl[int(100*tPS)+1000])
            APpPs = abs(rfl[int(100*tPpPs)+1000])
            APsPs_PpSs = abs(rfl[int(100*tPsPs_PpSs)+1000])
            
            A_sum += APS + APpPs + APsPs_PpSs
        Amp_sum[-1].append(A_sum)
xyzfile = '/Users/husir/Desktop/xyz.txt'
with open(xyzfile,'w') as f:
    for i in range(len(H)):
        for j in range(len(k)):
            f.write(str(H[i])+'    ')
            f.write(str(round(10*k[j])/10)+'    ') # 二进制存储导致有一些数字会产生偏移，这里是把一位小数后的数字都舍去
            f.write(str(Amp_sum[i][j]))
            f.write('\n')
# 存是存了个xyz文件，但是画图直接用了缓存着的几个数组，因为懒得重新读入还得调格式
import matplotlib.pyplot as plt
X,Y = np.meshgrid(H,k)
#print(X,'\n')
#print(Y,'\n')
Amp_sum=np.mat(Amp_sum).T
#print(Amp_sum)

# 绘制等值线图
plt.figure(dpi=200)
C = plt.contour(X,Y,Amp_sum,[1.2,1.5,2.0,2.5,2.7,2.9],colors='black')
plt.clabel(C,fontsize=6,colors='black')
#plt.contourf(X,Y,Amp_sum)
#plt.colorbar()
plt.xlabel('H/km')
plt.ylabel('Vp/Vs')
plt.savefig('等值线result')

最后结果是这样的：
等值线result.png

可以看到有一条斜下来的比较明显的条带上振幅值都比较高，最高的点是在H=30，k=1.8，虽然这个数值还算比较接近真实值，但是有这么多极大值点是我没有想到的，个中原因可能有大大小小的误差、各种操作的准确性等等。另外，这次计算是假定Vp=6.1km/s的，可以尝试多个别的Vp来对比最后的结果。

先写到这里，随时补充。
# 2020.7.11 19:42
btw 18分钟后Bilibili毕业歌会要直播了！我坐好了！

补充：
关于怎么从接收函数里拿到某个时间的振幅——道头信息里带有起始时间、终止时间和采样率，因此数据区域都是等时间间隔的振幅信息，只要将到时转换成数据数组的索引即可。
要注意的是，一般SAC文件以及我们转换出来的接受函数文件是二进制形式的，要能通过普通文本阅读器看到具体的道头信息、振幅信息，可以先用SAC转换成字符型形式。脚本如下：

#!/bin/bash
SAC_DISPLAY_COPYRIGHT=0
path=~/Desktop/data
cd $path
for eve in `ls`; do
    cd $path/$eve/wR_RF
    sac << EOF
      read ./decon.out
      write alpha rf  # 转换成字符型文件
      quit
EOF
    cat rf | sed '1,30d' > rf2.0         # 为删除前30行道头字符段
done

# 2020.7.12 09:30

补充：
saclst的用法也很简单，可以在终端中直接调用（没有想到要用的这么快）。saclst后面先跟要用的道头变量名，可以跟多个，列完了就用f表示后面是文件名了。

EVLO=`saclst EVLO f *BHE* | awk '{print $2}'`

# 2020.7.13 11:37

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 213,616评论 6赞 492
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 91,020评论 3赞 387
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 159,078评论 0赞 349
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 57,040评论 1赞 285
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 66,154评论 6赞 385
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 50,265评论 1赞 292
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,298评论 3赞 412
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 38,072评论 0赞 268
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,491评论 1赞 306
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,795评论 2赞 328
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,970评论 1赞 341
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,654评论 4赞 337
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 40,272评论 3赞 318
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,985评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,223评论 1赞 267
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,815评论 2赞 365
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,852评论 2赞 351