wangliang1989 / oh-my-cap

Oh My CAP

CAP 是计算震源机制的一种常用方法， 很多地震同行都需要学习、使用和研究它。 我基于 gCAP1.0 建立了 Oh My CAP 这个开源项目， 总结整理了在使用 CAP 中的经验， 希望能够为大家提供一个学习、探讨 CAP 的平台。 本项目始自 2016 年，我在桂林理工大学读硕士学位期间。

项目主页：http://seiswave.cn/oh-my-cap/

作者

王亮 贺州学院公共基础教学部

特别贡献者

田冬冬 **地质大学（武汉）

版权协议

本项目中所使用的 fk 和 gcap 的源码修改自朱露培的原始代码。 这部分代码遵循 GPL 协议，任何人均可免费获取、使用、修改和再发布代码， 但是修改后的版本必须也按 GPL 协议公开和授权。

本项目中的其余源码和文档采用更加宽松的 Apache 协议， 即在尊重本项目署名权的前提下，可以选择不公开自己的修改。

引引我的论文吧

如果本项目帮助了你，请引用我的论文。 有一些同学不愿意引用中文论文，或者不愿意引用非大牛的论文。 同学，讲点良心吧。 以下是我的论文的中英文信息：

王亮, 薛霆虓, 季海磊. 2016. 集集强余震震源机制解分析[J]. 地球物理学进展. 31(5): 1998-2004

Wang L, Xue T X, Ji H L. 2016. Focal mechanisms of Taiwan ChiChi earthquake aftershocks[J]. Progress in Geophysics. 31(5): 1998-2004

这篇论文可以在地球物理学进展的官网下载： http://www.progeophys.cn/CN/10.6038/pg20160515

另外，我也提供我的所有第一作者论文和学位论文。 在我的 坚果云网盘 中可以找到这篇论文，也可以找到我的以 CAP 为主题的硕士论文。

oh-my-cap's Issues

参数文件中参数名命名

数据预处理：
DT CAP里的-H参数，重采样的周期
FREQ 去仪器相应的四个频率

反演：
DEPTH 反演的深度
GREEN -G参数
MODEL 格林函数model
PLOT cap的-p参数，绘图
SHIFT cap -s
P_WINDOW 反演P窗口
S_WINDOW S窗口
WEIGHTFILE 权重文件
P_FILTER CAP -c参数里的滤波
S_FILTER 同上

以下三个我想不出合适名字：
-D
-W: 1
-X: 10

FK 计算修改配置文件规则方便修改模型

可能会更换模型，修改定义模型的方式，便于修改模型。另外需要确认给 fk 计算的模型文件一定会正确。

fk.pl检查fk文件的部分有错误，造成误报

foreach (@layer) {
        my $layer = trim($_);
        my ($dep) = split m/\s+/, $layer;
        foreach (@depth) {
            if ($_ == $dep) {
                print "设置错误：震源在 $_ km 深度的界面上\n";
                $err++;
            }
        }
     }

Mint安装

add rdseed

未考虑不同分量长度不同的情况

marktime.pl 未考虑不同分量长度不同的情况而造成 cut 数据出错

gmt6depth.pl报错

Ubuntu22.04 编译安装GMT6.1.1，运行示例perl gmt6depth.pl 20080418093658时报错
gmt [ERROR]: Shared GMT module not found: basemap

ERROR: No module named basemap was found. This could mean one of four cases:

1. There actually is no such module; please check your spelling.
2. You used a modern mode module name while running in GMT classic mode.
3. Module exists in the GMT supplemental library, but the library could not be found.
4. Module exists in a GMT custom library, but none was specified via GMT_CUSTOM_LIBS.
   Shared libraries must be in standard system paths or set via environmental parameter LD_LIBRARY_PATH.

gmt [ERROR]: Shared GMT module not found: plot

ERROR: No module named plot was found. This could mean one of four cases:

1. There actually is no such module; please check your spelling.
2. You used a modern mode module name while running in GMT classic mode.
3. Module exists in the GMT supplemental library, but the library could not be found.
4. Module exists in a GMT custom library, but none was specified via GMT_CUSTOM_LIBS.
   Shared libraries must be in standard system paths or set via environmental parameter LD_LIBRARY_PATH.

meca [ERROR]: Internal Failure = GMT_MAP_BAD_LAT_MAX
meca (GMT_psmeca): North is outside -90 to +90 degree range

需要完善的条目

1. 格林函数的计算，增加并行计算的脚本
2. 安装方法考虑为所用计算机账户安装的问题

将 hk 模型的格式修改为更常用的模型格式

目前 hk 中的四列数据分别为：厚度、vs、波速比、Qs，这不是一个常见的速度模型格式，建议修改为更通用的格式：

  thickness vs vp rho Qs Qp

下一个版本的新特性清单

1. 考虑到可能需要在工作区内放置一些文件，比如自己写的脚本，而又不想加入 repo 中，可
   以用 i_ 开头，则会忽略

gcap安装

老师您好，我是用里面的gcap时一直有问题
#报错：can‘t exec 'cap.pl':NO such...
请问是gcc要降到3.4版本的吗
另外，请问咱们是否有qq群这样的交流平台大家可以方便交流

如何降低波形拟合残差

请教老师，在导入自己的数据进行数据处理完之后，发现自己的波形拟合残差有点大，应该如何降低呢？

增加是否计算爆炸源格林函数的选项

现在脚本是按一律算爆炸源来处理，应该增加是否计算的选项

FK展平变换的说明

震源深度和模型速度的变换和地震学引论中一致。但是，模型界面深度和地震学引论中公式不同：
代码中的公式是 h'=r_0/(r_0-d+1/2h)，h是变换前的地层厚度，d是地层上界面的深度，r_0是地球半径，h'是变换后的地层厚度。
书中没有直接给出公式，可用深度的公式做推导：

权重脚本的修改

计划如下修改：
考虑加强权重设置的功能

gcap 外围脚本的重写

打算重写波形拟合图，拟合误差深度图和新写台站位置示意图

请教老师：运行例子的时候出现了一些问题！

增加批量处理事件的脚本

1. 新建脚本 build.pl: 读取事件目录文件，根据文件内容找到seed文件，然后建立对应的工作目录，把相应的事件信息保存为文本，将这个文本文件和seed文件移动到工作目录中去。
2. 新建脚本do_process.pl: 根据事件目录文件，去调用process.pl预处理各个事件。
3. 新建脚本do_inversion.pl:根据事件目录文件，去调用inversion.pl反演各个事件。

SAC 文件名、权重文件的协调问题

在反演时，是权重文件中的台名+分量来搜索sac文件的，对此要考虑有loc的时候，文件名协调的问题

确定例子事件

根据 gCAP 的官方版本中包含的例子中的sac文件的事件时间和经纬度信息，基本确认该事件是 IRIS 上的 2725882事件，但相关信息并不完全一致，请@seisman 确认

请教：gcap反演中不同时窗长度的P和S波，在输出图件中绘制的波形长度近乎一样的原因是为什么呢？

请教一下老师，设置不同时窗长度的P和S波，在输出图件中长度为什么几乎一样呢？在参数设置文件中，更改-P的值似乎也只能整体的改变其大小，如何才可以设置最终绘制的图件中各分量波形长度与其窗口长度是成比例的呢？谢谢。

参数的整理

1. 把 inversion.md 中对参数的介绍合并到 parameter.md 中，然后删除 inversion.md 中的部分，因为 inversion 的目的只是为了让读者对整个流程有一个认识而已。参数可以介绍，但没有必要这么复杂，想知道每个参数是什么含义就去看 parameters.
2. parameter.md 这个文件的格式要重新整理一下，现在的文档完全没有格式可言；
3. 搞不定的参数就空哪里空着，不要删除

请教：商业软件 Numerical Recipes

老师您好！请问gCAP需要使用的商业软件 Numerical Recipes在哪里可以购买

请教：去除仪器响应的问题

我在教程中看到：数据必须去除仪器响应得到位移或速度记录，单位为 cm 或 cm/s
我在使用seed文件中提取的PZ文件去除仪器响应后，（得到的结果据说单位是m）将结果乘以100，然后按照流程进行得到的反演结果看起来不太对。很多波形图超出边界，且计算出的震级与地震目录中的震级相差较大。
请问通过PZ文件去除仪器响应后是否有不一样的处理步骤？

关于 dt

共三个地方涉及到了 dt：

1. 格林函数
2. 数据
3. cap.pl 的 -H 选项

需要做的是：

1. 确认三个 dt 是否一定要相等
2. 若1成立，则修改源码使得其检测这一条件是否成立

配置参数化后更新文件的历史记录

配置参数化后更新文件的工作之前没有做完，做完后注意修改历史记录！

gCAP输出的结果

反演的矩张量结果保存在.out的结果文件中# MomentTensor那一行。
官方脚本是这样传入gmt绘图的：
$plt3 = "| psmeca -JX1/1 -R-1/1/-1/1 -Sm8 -G100 -Y9.5 -X-0.7 -O -K >> $outps" if $tensor[1] eq "MomentTensor";
open(PLT, $plt3);
printf PLT "0 0 0 @tensor[9,4,7,6] %f %f 17\n",-$tensor[8],-$tensor[5];
因为GMT的psmeca中，所用的是RTP/RTF/USE坐标系。所以得到对应关系（MOPAD默认用NED）：
9 4 7 6 -8 -5 = mrr, mtt, mff, mrt, mrf, mtf = dd nn ee dn -de -ne
或者说：
4 5 6 7 8 9 = nn ne dn ee de dd

请教：计算单力源和爆炸型震源时，程序报错。

请教一个问题：计算单力源和爆炸型震源时，程序报错。非常感谢！
命令行为：$ fk.pl -Mhk/15/k -N1024/0.1 -S0 -U0 10。（计算爆炸型震源时）。 程序默认的-S2（双力偶型震源）
报错：Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.
Backtrace for this error:
#0 0x7f8eac02532a
#1 0x7f8eac024503
#2 0x7f8eab8b9fcf
#3 0x559541734f9f
#4 0x559541734922
#5 0x5595417329db
#6 0x55954173380b
#7 0x7f8eab89cb96
#8 0x5595417313e9
#9 0xffffffffffffffff
hk_15/10.grn.0
Unable to open hk_15/10.grn.0
hk_15/10.grn.5
Unable to open hk_15/10.grn.5

删除不良数据

考虑在生成权重文件的脚本中加一段判断：如果 t9 不是 -12345，这个数据权重为0。

请教：输出结果中，波形图中无法出现P波的原因

运行过程中没有报错，但是在波形图中无法看到P波波形。
使用的事件波形的震中距在10~30km，震级均为2级多
请问老师，是什么缘故呢？

get_depth.pl 考虑分辨不同模型

get_depth.pl 调用 Zhu 的脚本，没有考虑到路径下可能有多个model的结果

需要确认psmeca命令作图存在bug，并寻找替代方案

#!/bin/bash
psmeca -P -JX0/10 -R-5/5/-5/5 -Sm8 -G100 -Y9.5 -X-0.7 > 1.ps <<EOF
0 0 0 -0.3 -1 -0.650 0 0 0 17
1 1 0 -0.3 -1 -0.649 0 0 0 17
EOF
ps2raster -E1080 -P -Tf 1.ps

使用以上内容作图，可以看到两个机制解其中一个黑色，然而这两个机制解应该均为白色。我这边作图结果是1.pdf

重写depth.pl并明确输出的含义

现在只实现了原depth.pl的做图功能，输出深度反演的功能还未实现

输出含义也未完全弄清

检查 gCAP 的做图脚本

1. 目前怀疑机制球台站位置绘制有错
2. 如果台站数很多，可能不能全部画上

以上问题需要检查

确认 TauP 是否必须

1. 确认 TauP 在整个流程中的作用
2. 若 TauP 非必须，则需要在 marktime.pl 中对是否安装了 TauP 进行检测。若未安装，则直接跳过标记理论到时