计算格林函数

使用 gCAP 反演震源机制的第一步是正确计算格林函数。 本项目默认将 Glib 作为格林函数库所在目录。

文件树如下：

Glib
├── model.fk                格林函数配置文件
├── config.pm               参数读取函数脚本
├── checkgrn.pl             格林函数波形绘制脚本
├── run_fk.pl               调用 fk，计算格林函数的脚本
├── run_parallel.pl         调用 fk，并行计算格林函数的脚本
└── model                   存放名为 model 的模型的格林函数的文件夹（在生成格林函数后才会有）
    ├── model               名为 model 的模型的文件
    ├── model_5             名为 model 的模型的震源深度为 5 km 的格林函数
    │   ├── [dist].grn.0    震中距为 dist 的 0 阶源产生的 Z 格林函数
    │   ├── [dist].grn.1    震中距为 dist 的 0 阶源产生的 R 格林函数
    │   ├── [dist].grn.2    震中距为 dist 的 0 阶源产生的 T 格林函数，振幅恒为 0
    │   ├── [dist].grn.3    震中距为 dist 的 1 阶源产生的 Z 格林函数
    │   ├── [dist].grn.4    震中距为 dist 的 1 阶源产生的 R 格林函数
    │   ├── [dist].grn.5    震中距为 dist 的 1 阶源产生的 T 格林函数
    │   ├── [dist].grn.6    震中距为 dist 的 2 阶源产生的 Z 格林函数
    │   ├── [dist].grn.7    震中距为 dist 的 2 阶源产生的 R 格林函数
    │   ├── [dist].grn.8    震中距为 dist 的 2 阶源产生的 T 格林函数
    │   ├── [dist].grn.a    震中距为 dist 的爆炸源产生的 Z 格林函数
    │   ├── [dist].grn.b    震中距为 dist 的爆炸源产生的 R 格林函数
    │   └── [dist].grn.c    震中距为 dist 的爆炸源产生的 T 格林函数
    ├── model_10            名为 model 的模型的震源深度为 10 km 的格林函数
    ├── model_15            名为 model 的模型的震源深度为 15 km 的格林函数
    ├── model_20            名为 model 的模型的震源深度为 20 km的格林函数
    └── model_25            名为 model 的模型的震源深度为 25 km的格林函数

格林函数配置文件

model.fk 是名为 model 的模型的格林函数配置文件。 这个文件的本质是一个文本文件，你可以用一般的文本编辑器打开。 fk 结尾的原因是这个文件用于 fk 计算格里函数。

model.fk 的内容如下：

# FK configuration
DIST: 135/145/5 210 235 260 280 300 415
DEPTH: 5 10/20/5 25
NT: 512
DT: 0.2
FLAT: NO               # NO or YES. if undefined, scripts will use YES

 5.5 3.18 5.50 2.53 600 1200
10.5 3.64 6.30 2.79 600 1200
16.0 3.87 6.70 2.91 600 1200
90.0 4.50 7.80 3.27 900 1800

说明如下：

  1. # 号后的内容对程序没有作用，只是注释；空行也没有作用
  2. DIST 是指格林函数的震中距，其中 135/145/5 等于计算 135、140 和 145 三个震中距
  3. DEPTH 是震源深度
  4. NT 是数据点数
  5. DT 是采样间隔
  6. FLAT 是设置是否进行展平变换
  7. 以 NT: 512 这一行来说，其中的英文冒号是必须的，冒号后的空格可以不要
  8. 程序支持的参数是 DIST、DEPTH、NT、DT 和 FLAT，没有其他的
  9. 最下面几行是一维地层模型，这部分的格式为：地层厚度、S波速度、P波速度、密度、S波Q值、P波Q值，后三项可以省略

如何正确设置计算的参数， 请参考 seisman 的 fk 用法笔记。

参数读取函数脚本

config.pm 包含了参数读取的若干函数。

计算示例格林函数

运行脚本 run_fk.pl 即可计算得到示例所需的格林函数。

$ perl run_fk.pl model.fk

如果安装了 Perl 的并行模块 Parallel::ForkManager， 你可以使用脚本 run_parallel.pl，用并行的方式完成上述计算：

$ perl run_parallel.pl model.fk

格林函数波形绘制脚本

checkgrn.pl 会读取格林函数的 SAC 文件，并将波形绘制为图片：

$ perl checkgrn.pl model.fk

pdf 格式的图片文件在路径 model 下。 在格林函数计算完成后，你应该检查格林函数的波形是否正确。 fk 已知的 bug 有当数据点数 NT 过短，导致后面的波形跑到前面去了。 所有，你需要看一眼。

格林函数库

不带后缀名的文件 model 是文本文件。 交给 fk 的就是这个文件。 计算完成后会得到 5 个名为 model_xx（xx 代表震源深度）的文件夹， 其中包含了震源深度取不同值所对应的格林函数。

以 model_5 目录为例，其中包含了众多文件名为 xxx.grn.x 的格林函数文件。其中：

  • xxx 代表震中距
  • x 取0-8，对应双力偶震源的9个格林函数分量
  • x 取a、b、c，对应爆炸源的3个格林函数分量

所有格林函数均是 SAC 格式， xxx.grn.0 和 xxx.grn.5 的头段 t1 和 t2 保存了初至 P、S 波的到时。 user1 和 user2 保存了初至 P、S 波的离源角（相对于垂直向下方向旋转的角度)。

gCAP 根据头段 t1 和 t2 的值确定要截取的波形时间窗； 若反演中需要使用震相初动极性信息，则需要从头段 user1 和 user2 中读取震相离源角。