简介

本文记录了利用 Oh-My-Cap 进行震源机制反演的基本流程。 按照教程进行安装,在安装 perl 模块的时候利用了 cpanm 进行安装。

拿到的数据是类似 AXI.SC.2017112153.00.BHE 这种命名格式的 SAC 文件,和 RESP.SC.AXI.00.BHE 命名格式的 仪器响应文件,把 SAC 数据文件夹和仪器响应文件夹放到 oh-my-cap/DATA 文件夹下。 建立自己的速度模型文件 qzgy.fk,执行脚本生成格林函数文件夹 qzgy。 根据要反演的地震信息修改 event.conf、event.info 和 RunMe.sh,运行 RunMe.sh 脚本即可开始反演。

# 在 DATA 文件夹运行 RunMe.sh 脚本进行反演
$ cd oh-my-cap/DATA
$ ls -1 
20170912.192642.000    # 地震数据文件夹
20170912.192653.000    # 地震数据文件夹
Catlog_cap.txt         # 地震目录
chosedata.sh           # 临时脚本,无用
draw_st.sh             # 绘制台站分布图
earth_relief_02m.grd   # 全球地形
event.conf             # cap 反演用到的配置文件
event.info             # 预处理用到的地震信息文件
gmt4depth.pl           # 王亮自己写的,还没用过不知道什么效果
GMT_globe.cpt          # cpt
mod.txt                 # 速度模型
my-process.pl           # 数据预处理脚本
rename.pl               # 修改命名格式
resp                   # RESP 仪器响应文件夹
RunMe.sh               # 反演的主脚本
sn.sh                # 计算信噪比
$ bash RunMe.sh

处理流程

计算格林函数

根据某地区的速度模型计算格林函数时可以把深度算的精细一点,然后根据需要在反演的时候选择适当的参数。

$ cat Glib/qzgy.fk
# 震源不能设置到分界面上,所以第一层层厚写成了5.1(或者4.9?)
# 只知道P波速度的情况下,Vs=Vp/sqrt(3)
# FK configuration
DIST: 15/500/5
DEPTH: 2/36/1
NT: 1024
DT: 0.2
FLAT: NO               # NO or YES. if undefined, scripts will use YES
# 空行内不能有空格,否则程序会将第一层的参数全部设置为0,导致无法正确生成格林函数
 5.1 3.23 5.6000
 5.0 3.34 5.8000
10.0 3.52 6.1000
12.0 3.63 6.3000
 8.0 3.81 6.6000
 6.0 3.92 6.8000
34.0 4.50 7.8000
$ cd Glib
$ perl run_parallel.pl qzgy.fk

event.info

查看地震目录文件或者利用 SAC 查看相关变量信息:

SAC> lh kzdate kztime b delta npts e o a f evla evlo evel evdp mag dist az

将相关信息填写到 event.info 中,填写的震级没有影响,但是深度对 Cap 有很小的间接影响, 因为 Taup 要依据深度算到时,算出来的到时用于 cut 波形,如果纯手工 cut 波形就没有影响。

# 这里需要注意时区问题,我拿到的地震目录中的发震时刻需要减去 8h
$ echo 2017-09-12T11:26:53 28.084 101.345 10 4.5 > event.info

event.conf

选择模型,会覆盖 example 文件夹下 event.conf 文件的配置信息。

$ cat event.conf
# CAP paramete# CAP parameters
MODEL: qzgy
DEPTH: 3/36/3
MAG: 4.0

# CAP options
-J: 0/0.1/0/0.1

数据预处理

执行 process.pl 报错,查看脚本源代码发现需要按照固定的格式命名文件,编写 rename.pl 脚本进行重命名

#!/usr/bin/env perl
use strict;
use warnings;
$ENV{SAC_DISPLAY_COPYRIGHT} = 0;

@ARGV == 1 or die "Usage: perl $0 dir\n";
my ($dir) = @ARGV;

chdir $dir;

# 文件重命名
foreach my $file (glob "*20*.BH?") {
    my ($sta, $net, $time, $loc, $chn) = split /\./, $file;
    rename $file, "$net.$sta.$loc.$chn.SAC";
}

chdir "..";

将 process.pl 改成 my-process.pl (其实就是把里面的 rdseed.pl 替换为 rename.pl),执行脚本

$ perl my-process.pl chosedata

反演

# 必须先生成权重文件
$ perl weight.pl 20170912.192653.000
# 反演
$ perl inversion.pl 20170912.192653.000
# 深度反演
$ perl get_depth.pl 20170912.192653.000

如果出现类似 Warning: flag=21 => the minimum 0.0/57.8/ 6.4 is at boundary 的警告,证明某个台站位于节面上, 特点是 P 波非常的弱,在权重文件中去除这个台站就可以。当然,也可以选择忽略这个警告,对结果没有太大影响。

如果出现类似 no minimum found 的警告,可能原因有:

  1. 某个台站的波形经过预处理后“消失了”(目前没查到原因),利用 SAC 绘图找到这个台站后在权重文件中删除这个台站即可。目前发现利用 RESP 去除仪器响应波形数据消失的情况下改用 PZ 就不会出现问题。

  2. 格林函数没有正确生成。请注意qzgy.fk参数配置文件中速度结构部分空行内不要有空格。

绘制台站分布图

读取 event.info 中的震中经纬度,绘制经纬度 ±5° 范围内的台站分布图。

#!/bin/bash
awk '{print $3,$2}' event.info|awk '{
                lon1=int($1-5);lon2=int($1+5);
                lat1=int($2-5);lat2=int($2+5);
                print "grdcut *.grd -R"lon1"/"lon2"/"lat1"/"lat2" -G1.grd";              
                print "grdgradient 1.grd -Gint.grd -Nt1 -A30/60";
                print "grdimage -R -JM5i -CGMT_globe.cpt -Iint.grd 1.grd -K -P -B1 > st.ps";
        print "psxy -Sa0.3i -R -J -Gred -K -O >> st.ps <<EOF";
                print $1,$2;
                print "EOF";}'|sh
pscoast -R -J -W1p -K -O >> st.ps
awk '{print $2,$3}' st_info|psxy -St0.2i -R -J -G55/0/250 -K -O >> st.ps
awk '{print $2,$3+0.1,"12 0 0 BC",$6}' st_info|pstext -R -J -O >> st.ps
ps2raster st.ps -A -Tf
rm .gmt* *.grd st.ps 

计算信噪比

利用黑板变量和头段变量计算每个台站 Z 分量的信噪比。 采用的信号范围为 Pnl 波到时(t0)之后的 100s,采用的噪声范围为数据结束(e)之前的 100s。 信号没有进行噪声校正。

#!/bin/bash
echo \# stlo stla az dist s/n kstnm > sn_info
ls -1 *.z|awk '{print "r",$1;
                print "setbb var1 (&1,t0&)";
                print "setbb var2 (&1,t0&+100)";
                print "setbb var3 (&1,e&-100)";
                print "setbb var4 (&1,e&)";
                print "mtw %var1% %var2%";
                print "rms to user5";
                print "mtw %var3% %var4%"
                print "rms to user6";
                print "setbb sn (&1,user5&/&1,user6&)";
                print "ch user7 %sn%";
                print "lh stlo stla az dist user7 kstnm";
}END{print "q";}'|sac| \
           awk '/ stlo /{printf "%-8.4f ",$3}
        / stla /{printf "%-8.4f ",$3}
        / az /{printf "%-8.4f ",$3}
        / dist /{printf "%-8.4f ",$3}
                / user7 /{printf "%-8.4f ",$3}
        / kstnm/{printf "%s\n",$3}' | sort -k 5 -rn >> sn_info

说明

  1. Cap 研究中需要的频段是 0.05 到 0.2 Hz,f1|f2|f3|f4 可以取为 0.02|0.03|0.5|1,这个基本是最小的范围,也可以取更大范围,比如 0.01|0.02|20|30,因为 trans 是第一次滤波,cap.pl 中 -C 选项进行了第二次实际滤波,另外,地震的能量主要部分也不在高频,震级越小,信号越高频,所以小震机制解是难点。

  2. 初始速度模型很关键,台站的选取也影响到最后的结果。可以选取四象限分布较好的震中距范围在 100-300km 范围内的几个信噪比较好的台站进行反演,因为震中距太大波形不清晰,震中距太小面波没到。

  3. 反演结果中 ISO、CLVD 代表非构造地震的成分,默认不反演,可以在 event.conf 中用 -J: 0/0.1/0/0.1 控制。 深度反演结果中横坐标为深度,纵坐标为误差,h 后的两个数为拟合结果最好的深度以及深度估计范围。

附录-RunMe.sh

#!/bin/bash
begin_time=$(date +%s)
# 设置需要反演的目标文件夹(需要手工修改)
dir=20140105.202145.000

if [ ! -d $dir ]; then
  echo ===== 文件夹 $dir 不存在,终止程序! 
  exit 1
fi

echo ===== 目标文件夹为 $dir
rm -rf ../$dir
# 根据地震目录修改两个文件信息后复制到目标文件夹(需要手工修改)
echo ===== 正在复制 event.info 和 event.conf 到 $dir
cp event.conf event.info $dir
echo ===== 正在复制仪器响应文件到 $dir
cp resp/* $dir
cd $dir
rm -rf chosedata
mkdir chosedata
echo ===== 正在使用 saclst 提取震中距 80-300km 的台站,按照方位角降序排列
saclst stlo stla az dist kstnm f *20*.BHZ|awk '$5>80 && $5<300 {print $0}'|sort -k 4 -rn > st_info
echo ===== 正在将筛选出的台站和台站的仪器响应文件等复制到 chosedata 文件夹
awk -F . '{print "cp event.info event.conf st_info ""*"$1"* chosedata"}' st_info|sh
cd ..
echo ===== 正在将 chosedata 文件夹复制到 example 文件下并改为 $dir
cp -r $dir/chosedata ../$dir
echo ===== 正在进行数据预处理
cp my-process.pl rename.pl ../
cd ..
perl my-process.pl $dir
echo ===== 正在生成权重文件
perl weight.pl $dir
echo ===== 正在进行反演
perl inversion.pl $dir
echo ===== 正在进行深度反演
perl get_depth.pl $dir
rm my-process.pl rename.pl
echo ===== 正在绘制台站分布图
cp DATA/draw_st.sh ./$dir
data=DATA/earth_relief_02m.grd
cpt=DATA/GMT_globe.cpt
cp $data $cpt ./$dir
cp DATA/sn.sh ./$dir
cd $dir
bash draw_st.sh
echo ===== 正在产出台站信噪比文件 sn_info
bash sn.sh
end_time=$(date +%s)
declare -i cost_time=$end_time-$begin_time
echo ===== Done! =====
echo ===== 本次反演总共花费 $cost_time 秒 =====

参考资料