[Literature Reading]Pytheas: An open-source software solution for local shear-wave splitting studies

[Lszidv]

Abstract

提供了包括视觉检查、旋转相关法、特征值法和最小能量法在内的多种分析工具,并通过聚类分析实现自动化处理。
支持标准数据格式，如QuakeML、StationXML、SAC和miniSEED
引入了一种新算法，用于自动为测量结果分配质量等级
手动分析建立对照组，完全自动化处理数据目录，并通过视觉检查样本来验证结果。

1. Introduction

剪切波分裂（SwS）介绍及其在地壳各向异性研究中的重要性。

• 事件台站配对数据集，数据被拒绝主要是因为剪切波难以识别（噪声、次生波相的干扰）以及次生波相（转换波或散射波）的干扰
  为了减少次生波的影响，限制仅考虑在剪切波窗口内传播的射线（利用区域的速度结构来设定入射角度）
  ES系统导致一些数据被误筛

在粒子运动图中识别ϕ量，对波形进行各向异性矫正后，寻找最好的时间延迟δt
轨迹图（hodogram）, 极化图（polarigram）

原始水平分量（E-W和N-S)
- 原始水平分量是地震仪记录下的地震波的水平分量
- 这些原始分量是地震数据的基础，它们为后续的分析提供了波的振幅、频率和传播方向等基本信息。通过原始水平分量，可以进一步推导出径向、切向分量，以及快慢波分量。
---
径向和切向分量（Radial and Transverse Components）
- 径向分量：指沿着震中和观测点连线方向的波动分量，即从震中传播到观测站的方向。
- 切向分量：指垂直于径向分量方向的波动分量，即水平面内与震中连线垂直的方向。
- 径向分量主要用于捕捉波沿震中到台站方向的传播特性，而切向分量则用于描述波在水平方向上垂直于震中连线的传播情况。
-径向和切向分量之间的转换通过坐标旋转实现，旋转的角度通常基于震中和台站之间的几何关系。
---
Sfast和Sslow分量（Sfast and Sslow Components）
- Sfast（快波）：是剪切波分裂现象中的一个分量，传播速度较快，沿着Sfast的偏振方向传播。
- Sslow（慢波）：是剪切波分裂现象中的另一个分量，传播速度较慢，沿着Sslow的偏振方向传播。
- 通过对径向和切向分量进行坐标旋转（ϕ）来获得（原始分量和快慢波的偏振方向不一致，R、T分量作为中间步骤，提供了更好的基准方向）。旋转后的坐标系统使得Sfast和Sslow分量能够沿其偏振方向进行区分，从而分析它们的传播速度差异。

网格搜索法：最小能量法（网格搜索寻找能最小化切向分量能量的最佳参数对）和特征值法（评估径向和切向分量的协方差矩阵的第二特征值，最小化这个特征值。）
旋转相关法：将水平波形旋转到不同的F-S坐标系中，并为每个旋转后的坐标系应用不同的时间延迟值（td）。
特征值法（EV）与聚类分析（CA）相结合：
  1.使用EV方法，首先为围绕Sfast到达的信号窗口获取多个参数对（ϕ和td）；
  2.根据参数对之间的接近度，在ϕ和td的二维空间中进行聚类分析；
  3.最终，通过选择最“受约束的聚类”来确定最稳定的解，从而获得最优的分析结果。

pytheas对以上方法进行集成，整合到一个图形用户界面中，并且内置了数据库

2. The Pytheas software

使用Obspy（用于地震数据处理）和Matplotlib（用于生成图形）。Scikit-learn（用于聚类分析）
Qt5框架和PyQT5绑定确保了Pytheas拥有易于使用且跨平台的图形用户界面。
使用SQLite3来存储分析结果，避免了复杂的服务器设置，支持本地存储。
2.1. Software design
一个脚本六个模块

主脚本（pytheas.py）：负责构建图形用户界面（GUI）并将各种功能连接到界面上。
六个模块：
parsers：用于读取和写入配置文件的设置。
db_handler：处理与数据库的所有交互。
eigenvalue：实现特征值法（EV）分析。[特征值法（EV）和最小能量法（ME）所需的函数（Silver和Chan，1991）;ϕ和δt的误差公式，Walsh等（2013]
rotationcorrelation：实现旋转相关法（RC）分析。[Bowman和Ando（1987）；通过计算95%置信区间，Pytheas为参数（如ϕ和δt）提供了一个误差范围。]

数据库包含三个表：
event：存储事件元数据。
station：存储台站元数据。
method：存储与SwS测量相关的信息，如ϕ和td值。

clustering：执行聚类分析（CA）。
Teanby等（2004b）提出的聚类分析（CA）方法
聚类模块使用一个类对象来管理分析流程中的各个步骤：
   候选窗口生成：定义用于分析的信号窗口范围
   聚类形成：根据相似性对数据进行分组
   聚类选择：选择最优的聚类结果，用于进一步分析

左侧窗口起始点：从S-P时间（S波到达台站的时间减去P波到达台站的时间）的一半开始，这有助于避免选择过短的窗口。※[Savage等（2010）]
右侧窗口结束点：由S波到达时的主频率 * 自定义的因子(对右侧结束点的时间进行扩展或缩小)

层次聚类：递归地将数据点合并为更大的聚类，或将大的聚类分割为小的聚类，最终形成一个树状结构
四种链接聚类的标准，都使用欧几里得距离（直线距离）作为相似性度量：
ward标准：最小化合并后的聚类的方差。
average标准：最小化两个聚类之间的平均距离。
complete标准：最小化两个聚类之间的最大距离。
single标准：使用最小的距离来合并聚类。

tools：提供辅助工具和功能。
自动评分算法：
  1.检查数据是否有效：检查ϕ与修正剪切波的偏振方向的角度差（无效记为 ‘N’ ）
      定义无效测量：没有剪切波分裂现象或S波的偏振方向过于接近快波或慢波的偏振方向（自定义阈值标准，用于判断近似平行与垂直）
  2.使用最小信噪比筛掉有显著噪声的记录（‘E’）
  3.评分：Pytheas的自动评分算法通过计算四个单独的评分，来评估剪切波分裂分析中的测量结果质量。具体评分指标包括：
    ϕ的误差（S-δϕ）；
    时间延迟（td）的误差（S-δt）；
    FS轴系统中的波形相关性（S-CC-FS）；
    NE轴系统中的波形相关性（S-CC-NE）。
每个评分的计算方式是基于实际测量值与预定义的最大和最小阈值之间的差异。通过这些评分，系统可以判断测量结果的可靠性，并根据评分值判断数据是否有效。如果某个评分值超过1.0，测量会被拒绝（标记为“E”级）。如果所有评分都在允许范围内，则计算总评分S，并通过平均值来决定数据的质量。
    评分由D到A,每个评分参数的阈值和等级都可自定义，为计算误差，必须自动评分。
额外功能：
  Obspy的TauP：该工具用于计算理论上的S波到达时间和入射角。TauP方法用于模拟波的传播路径，并根据给定的介质结构计算S波到达观测台站的时间和入射角。
  AR-AIC自动拾取器：这是一个自动化工具，用于从地震记录中自动提取S波到达时间。**AIC（赤池信息准则）**是一个常用的信号处理方法，用于选择最适合的模型或参数。

（1）. sδϕ（ϕ的误差评分）：

$$ s_{\delta \phi} = \frac{\delta \phi}{\delta \phi_{\text{max}}} $$

• 其中，δϕ 是偏振方向ϕ的误差，δϕmax 是允许的最大误差。

（2）. sδt（时间延迟误差评分）：

$$ s_{\delta t} = \frac{\delta td}{\delta td_{\text{max}}} $$

• 其中，δtd 是时间延迟（td）的误差，δtdmax 是允许的最大时间延迟误差。

（3）. sCCFS（FS轴系统中的相关系数评分）：

$$ s_{\text{CCFS}} = \frac{1 - |CCFS|}{1 - CCFS_{\text{min}}} $$

• 其中，CCFS 是FS轴系统中修正后的水平波形的相关系数，CCFSmin 是允许的最小相关系数。

（4）. sCCNE（NE轴系统中的相关系数评分）：

$$ s_{\text{CCNE}} = \frac{1 - |CCNE|}{1 - CCNE_{\text{min}}} $$

• 其中，CCNE 是NE轴系统中修正后的水平波形的相关系数，CCNEmin 是允许的最小相关系数。

2.2. Analysis workflow
软件三个必须输入：
事件目录文件：包含地震事件的详细信息；[QuakeML文件或简单的ASCII文件]
台站信息文件：包含台站的相关信息波形数据存储路径；[StationXML或纯文本文件]
指向存储地震波形数据的文件夹路径。[波形文件：Obspy可读取的任何格式（例如，SAC和miniSEED）]
（每个事件记录存储在以震源时间命名的文件夹中，每个台站记录存储在单独文件夹中）
结合FDSN（国际数字地震仪网络联合会）的标准化接口

2.2.1. Manual method

挑选Sfast波到达时间：分析师手动确定Sfast波的到达时间。
获取极化方向：软件从NE平面的运动向量中获取Sfast波的极化方向ϕ。
旋转到FS轴向系统：将水平分量旋转到FS轴向系统。
确定时间延迟：通过时间上移动Sslow波，以匹配两个分裂波的到达时间，从而确定时间延迟td。
旋转回NE系统：将波形旋转回原始的NE系统。
自动确定极化方向：Pytheas软件自动确定校正后的剪切波的极化方向。
手动提供等级：分析师需要手动提供一个等级，以评估结果的质量。
生成总结图表：可以生成一个图来总结分析结果。

2.2.2. Semi-automatic methods
选择分析窗口（RC、EV 和 ME）

左上：初始NE坐标系中经过滤波后的水平波形；旋转到FS坐标系并经过各项异性矫正；旋转回NE坐标系并再次各项异性矫正
左下：初始以及经过经过矫正的粒子运动图
右侧：特征值轮廓图，黑色区域为95%置信区间

2.2.3. Automatic method
使用CA方法，有两种执行方式 ：
1.用户选择台站事件配对，手动或者从事件记录中提取Sfast到达时间，然后选择方法组合进行分析(如CA-RC、CA-EV或CA-ME)
2.自动分析数据集中所有可用的台站-事件配对，但到达时间必须在QuakeML文件中提供（可以选择剪切波窗口）
生成聚类质量评估图与上述半自动方法质量评估图

初始数据的分布：展示了候选信号窗口的位置和分布。
Calinski-Harabasz评分：展示了不同聚类数目的质量评估，帮助选择最佳聚类数目。
聚类结果：展示了数据的聚类结果，将相似的信号窗口分到同一组。
聚类和窗口选择：展示了选择最合适聚类和信号窗口的过程。
ϕ变化：展示了不同窗口下Sfast波偏振方向的变化。
td变化：展示了不同窗口下Sfast波和Sslow波之间时间延迟的变化。

3. Case study: the Western Gulf of Corinth
   Western Gulf of Corinth (WGoC)一个月的地震记录，HUSN和CRLN两个台阵
   263个地震事件，震级范围从0.6到4.0，分析包括了一个Mw 5.1的地震事件
   (Papadimitriou et al., 1999; Kaviris et al., 2017, 2018b)，剪切波窗口被设置为45°
   CA-EV方法（偏好设置见表）进行自动分析，并且通过MAN方法分析了56个事件作为对照组。
   所有数据的相位到达时间和事件目录来自雅典大学地震实验室
   

评估过程：
确定了136个无效测量结果，其中“E”等级的结果为930个，“D”等级为20个，“C”等级为146个，“B”等级为177个，“A”等级为30个
无效测量有一些是台站N-S分量未工作，“E”等级的结果质量较差，其中有158个测量的时间延迟（td）过高，超过200毫秒
比较两种方法共有结果，仅接受C-A ，使用皮尔逊相关系数（PCC）

φ---0.71 δt---0.53

4. Conclusions

结合现有分析方法，自动化，结果可视化
提供与FDSN的接口
集成的数据库便于生成玫瑰图和φ与δt的时间变化图

未来方向：
远震
加权评分机制
滤波选择和优化 Savage et al., 2010
并行计算

## 自动化功能：
1. **数据预处理**：自动下载、检查、切割地震数据并计算信噪比。
2. **时间窗口选择**：基于波形特征自动选择合适的时间窗口。
3. **剪切波分裂分析**：自动化地计算分裂时间和极化方向。
4. **模型拟合**：自动拟合单层或双层各向异性模型。
5. **误差计算与优化**：自动计算误差，并优化结果的稳定性。

## 手动干预：
1. **质量检查与用户确认**：用户需要手动检查和修改自动选择的时间窗口。
2. **数据剔除与干预**：需要手动确认是否将某些空测量数据包含在分析中。
3. **结果评估**：用户需要评估和确认分裂结果的质量，尤其是在复杂波形和噪声较高的情况下。

Appendix

偏振方向玫瑰花图对比示例

Paper Link