class RFStation

[JouSF-1409]

最近有空琢磨源码的组织，想提升下泛用性，起因是当下使用的数据确实有些复杂，也不太想动文件的组织形式，数据量较大，担心改名等操作的时候出些问题。

最主要的希望是文件管理上能有更大的自由。
比如说传入一个命名标准，在调用时由sta_lst和其他的信息生成文件名。

读文件、处理数据的几个功能做一定的分离，改进一下接口（这部分主要是在ccp3d这里）

然后是S波接收函数。正好在做，所以像拓展一下ccp3d的适用性。

[JouSF-1409]

我还没读过代码的其他部分，如果我参照obspy.core.trace.Trace的感觉去实现，在ccp3d或者其他脚本中是否会有很多问题。

[xumi1993]

最近有空琢磨源码的组织，想提升下泛用性，起因是当下使用的数据确实有些复杂，也不太想动文件的组织形式，数据量较大，担心改名等操作的时候出些问题。

最主要的希望是文件管理上能有更大的自由。 比如说传入一个命名标准，在调用时由sta_lst和其他的信息生成文件名。

读文件、处理数据的几个功能做一定的分离，改进一下接口（这部分主要是在ccp3d这里）

然后是S波接收函数。正好在做，所以像拓展一下ccp3d的适用性。

• ccp3d只需要一个输入文件（时深转换之后的数据体ccppara.CCPPara.depthdat），原则上不存在文件管理方面的问题。方便举个例子详细说明吗？
• S波接收函数我没怎么测试过，可能时深转换的部分也需要debug，如果能测试一下就十分感谢了🙏

[JouSF-1409]

时深转换之前的步骤是用Prof. Zhu 的iter_decon做的。一部分文件按 事件/文件的形式整理，另一部分刚从seed解压出来。所以文件名命名有点混乱。
之后的步骤没能做下去，一部分是因为算穿刺点的时候，旧的脚本出现了各种各样的问题。另一部分是因为做LAB并不是410-610 or Moho，深度变化比较大。想引出bootstrap叠加筛选gui的那个接口看看怎么用。

所以我算是从中间开始，并没有严格按照seispy的习惯来做。所以很多东西需要人为去对齐，包括
文件名，射线参数的单位等。

rf2depth_makedata.py和rfcorrect.py这里前后端没有很明确地分离(相比之下，obspy.core.stream的实现要干净一些，可能我脑子还是在用C去理解python)。这里面stadatar，sta_info，两个变量，class Station，class RFStation几个在用法/含义上有些重叠。
重写的时候会按我这边更一般性的要求，参考obspy源码（引入doctest，processing_info等）写一些功能，不过对其他用到RFStation的地方，像是H-k或是ccp_profile等地方就暂时雇不到了。

在做S波接收函数之前，想先了解的几件事：

1. 静校正，ustclitho2.0除了常规的模型，还提供了一个以水准面为标准的速度模型，其高程最高可到-7。seispy能做静校正，还是最后出结果的时候叠加地形图再看。
2. 使用三维速度模型计算S波接收函数穿刺点和时深转换要对速度模型或者参数做哪些调整。

[xumi1993]

• Prof. Zhu 的程序CCP叠加是一步，而seispy是rf2depth和ccp_profile/ccp3d，所以如果要用iter_decon计算的RF，必须同时做时深转换和CCP叠加，可以写个脚本转换成seispy的格式，这一步可以单独提供一个脚本，在文档中说明。
• rfcorrect.py里是class，所有单一台站的RF处理都调用了这个class，rf2depth_makedata.py是函数也调用了rfcorrect.RFStation。rf2depth_makedata.Station只是用于读取台站列表，没有其他功能，可能在命名上有歧义。
• 我不建议先修改rfcorrect.py，因为他可能导致其他功能的损坏，可以先修改rf2depth_makedata，把makedata写成class。
• doctest是一个非常好的模块，但我担心seispy很多函数功能比较复杂，有的涉及数据，实现比较复杂。seispy很多函数缺少docstring，这是一个长远的任务。
• 我不清楚processing_info是什么模块。

[xumi1993]

在做S波接收函数之前，想先了解的几件事：

1. 静校正，ustclitho2.0除了常规的模型，还提供了一个以水准面为标准的速度模型，其高程最高可到-7。seispy能做静校正，还是最后出结果的时候叠加地形图再看。
2. 使用三维速度模型计算S波接收函数穿刺点和时深转换要对速度模型或者参数做哪些调整。

理论上S波的时差校正和时深转换和P波是一样的，seispy可以处理带高程的1D速度结构，最后按照给定的深度范围插值，可以参考这个class

seispy/seispy/utils.py

Line 70 in 31b67e3

class DepModel(object):

[JouSF-1409]

• Prof. Zhu 的程序CCP叠加是一步，而seispy是rf2depth和ccp_profile/ccp3d，所以如果要用iter_decon计算的RF，必须同时做时深转换和CCP叠加，可以写个脚本转换成seispy的格式，这一步可以单独提供一个脚本，在文档中说明。

您说的这些内容，我确实读到了。目前的计划是自己建类，然后再看情况拓展。
主要的目的是把io剥离开类的实现，原因有2：

1. 文件类型，命名方式，参考时间设置，不同来源的不太一样。在这一方面似乎应该留出足够的自由度。类似于obspy.io的设计。
2. python的glob等方法，在执行效率上似乎有一定问题。之前写过一个用glob遍历，obspy统计头部信息的脚本。但运行一段时间后运行速度明显下降，可能与操作系统的调度有关。使用一个文本文档做 索引能显著改善这个问题。这也是我想做的东西。

doctest，在我读的obspy源码中，还只是在处理字符串，比如说UTCDateTime类的实现中用的比较多。obspy中类似的工具与设计还有很多，注释也写得很完整，大半的本领都来自于此。

processing_info是一个decorator，在使用obspy.core.trace.Trace对象的内置方法(如trim)时，能自动向Trace.processing中添加调用时的func.__name__与参数，类似于shell中使用上方向键。向seispy中引入的目标是在运行时向setuplog中传入必要信息。

DepModel是一个很好的功能，尤其是现行的不同软件，对速度模型文档的结构有不同的要求。这样一个类减轻了很多负担（尽管我还没细看其实现）。

我不是很清楚您说的内容，可否理解为，倒转的S波接收函数序列，从射线的角度说，可以直接套用P波接收函数的计算方法(自己的琢磨。

[xumi1993]

主要的目的是把io剥离开类的实现，原因有2：

1. 文件类型，命名方式，参考时间设置，不同来源的不太一样。在这一方面似乎应该留出足够的自由度。类似于obspy.io的设计。

• 对于这部分，我并不希望对seispy添加这样的功能，seispy和obspy有一个本质区别是，seispy侧重于数据处理流程，obspy主要提供API，我希望能保证从接收函数计算到后续处理的数据的规范性，我不希望在中间步骤扩展其兼容性。与obspy不同，这一步只是文件存放方式的差别，而不是定义了另一种文件格式。

2. python的glob等方法，在执行效率上似乎有一定问题。之前写过一个用glob遍历，obspy统计头部信息的脚本。但运行一段时间后运行速度明显下降，可能与操作系统的调度有关。使用一个文本文档做 索引能显著改善这个问题。这也是我想做的东西。

• 如果有好的方法，可以改进，但我用文本文档做索引显然不是pythonic 习惯。而且目前在这部分我没有遇到glob执行效率的问题，或者说目前的效率我是可以接受的，我最多计算过4000+接收函数没有感觉到明显的效率问题。CCP叠加主要慢在搜索每个bin里的穿透点和bootstrap，这一步如果可以用JIT语言，应该能提高很多。
• doctest是用docstring中的命令做集成测试用的，避免了文档和测试的重复开发，目前seispy的集成测试脚本是单独写的 seispy/test。 这些测试都用到了数据，所以doctest可能不适用所有的函数，但是仍然有必要把所有函数的docstring补全。

processing_info是一个decorator，在使用obspy.core.trace.Trace对象的内置方法(如trim)时，能自动向Trace.processing中添加调用时的func.__name__与参数，类似于shell中使用上方向键。向seispy中引入的目标是在运行时向setuplog中传入必要信息。

• 这确实是一个不错的想法，你可以针对这个问题提一个issue，逐步推进这个功能。

[JouSF-1409]

我又回来了
这次讨论下 DepModel里的radius_s。计算穿刺点piercing point，或者说偏移位置的方法。

计算采用的公式

做了一些弧长与高的近似。但这个三角形，我是指圆心-上下两层的交点形成的三角形，是能用正弦定理表示并解出每个夹角的准确值，不需要做进一步近似。

这里是有什么特别的考虑吗。

[xumi1993]

高是指什么，$i$ 是第 $i$层不是角度 ，那 $coti$是什么。我猜你的意思是每层的入射角不一样，可以算出来。这样计算就是对射线路径的校正。 随着射线的偏移，速度横向会产生变化，也要考虑进去的话，每条射线的每个网格都要算一次角度，这个计算量是非常大的，用r f2depth -r可以实现这个功能，但计算量相当大，而且这依然是基于层状模型的，入射角也不准确。射线不确定性的影响也要远远小于速度的影响。

[JouSF-1409]

凌晨的时候不太清醒。抱歉。之前有文章讨论过平面波假设的问题，影响约在3%以下（660界面）。那么对于大部分界面来说，转换波与直达波之间射线参数的变化不应当为误差的主要因数。

简单画一下图，theta和j表示角度。上面的简化是将p=R*sin j /v 展开后得到的，那么这个近似就比较直观了。就是这个垂线的长度与R的比作为角度，这个近似可以几种不同的说法，anyways

但如果说用三角公式去算，有

既知 j, r1, dz 等，那么就能用反三角函数去算 theta。考虑到 sin在 0~pi的对称性，可以转换到cos再计算。

[xumi1993]

这样也可以，但影响应该不大。

[JouSF-1409]

我有点迷惑comp在匹配文件 finallist.dat还有相关流程中的作用。

像是rfcorrect.py里的 RFStation里定义了 only_r， read_stream也需要一个prime_comp。
如果只是跑一边常规的流程，comp似乎都不需要定义。似乎是研究各向异性时会用到？

然后关于之前说的add_processing_info这个修饰器。实现是比较简单的。但是 在哪些函数上添加是个问题，确实能返回参数与值，或者输出到日志系统，像下面这样。但是如果调用的是numpy.ndarray或者其他结构比较复杂的数据就比较复杂了。

日志的输出上，还是选择使用类变量做的。我觉得这样做比较节省创建日志的开销？ 毕竟DepMod和RFStation还是会经常调用到的。

[JouSF-1409]

我还有一件比较奇怪的事情是，RFStation 的data_prime使用二维数组存储波形数据，
但是我看到的脚本好像都没用到二维数组的优点？都是用ev_num做range进行迭代

[xumi1993]

comp是留下了的一个接口，因为理论上用PP PKP也可以算接收函数。
你可以把每一个要添加的功能先提交一个issue，进行针对性讨论。

[xumi1993]

我还有一件比较奇怪的事情是，RFStation 的data_prime使用二维数组存储波形数据，

但是我看到的脚本好像都没用到二维数组的优点？都是用ev_num做range进行迭代

这是因为这些功能都不能直接进行矩阵运算。