Abstract

1.Introduction

1. SWSPy简介

SWSPy专门用于剪切波分裂分析，能够准确高效地测量单个地震-接收器射线路径的剪切波速度各向异性。它使用Python编写，易于安装和集成，支持并行化处理。

2. 自动化与大规模数据处理

SWSPy设计为完全自动化，能够处理数千个地震事件和接收器的大规模地震数据集。随着地震仪器和数据存储技术的进步，SWSPy可以处理更大规模的数据集，从而增加检测到的地震数量。尽管SWSPy使用Python编写，但其计算密集型部分采用编译语言实现，以提高效率，优化软件性能。

3. 支持三维与多层分裂测量

SWSPy支持三维剪切波分裂测量（Walsh等)，并能在某些情况下分析多层剪切波分裂测量，特别适用于复杂地质条件下的研究。

4. 与其他软件的对比

SWSPy与现有的其他剪切波分裂分析软件相比具有显著优势。许多现有的软件包，如MFAST（Savage等，2010）和SHEBA（Wuestefeld等，2010），通常要求用户手动输入大量数据，并手动进行数据窗的选择。尽管这些软件包已经在学术界得到广泛应用，但它们通常缺乏自动化处理的能力。SWSPy则通过完全自动化的设计，在处理大规模数据集时大大简化了操作流程。此外，SWSPy的并行化设计可以充分利用现代计算机硬件的优势，进一步提升数据处理能力。与SplitLab、Pytheas等软件相比，SWSPy支持多层剪切波分裂分析，能够更全面地处理复杂的地震数据。

5. 软件性能展示

在本研究中，展示了SWSPy的多种应用实例，证明了该软件在大规模数据处理中的实用性和性能。

2.methods

在各向异性介质中的剪切波分裂可以通过两个主要参数来描述：

为了测量这些参数，存在多种方法，包括：

  1. 互相关法：通过旋转和平移射线垂直（通常为水平）分量，最大化两者的互相关，来找到最佳的分裂参数。
  2. 分裂强度法：通过接收器处多个地震事件的横向分量的本征向量的方位依赖性来确定分裂参数。该方法要求延迟时间小于剪切波的主周期，并且需要有足够的后方位角覆盖。
  3. 特征值法：通过旋转和平移射线垂直分量，寻找最小的特征值比率，从而确定剪切波分裂参数。该方法在大多数剪切波分裂应用中稳定可靠，并且适用于各种局部到远震的数据。

SWSPy软件包采用了特征值法，并结合了 Teanby et al. (2004) 的多窗口聚类技术和 Walsh et al. (2013) 定义的三维坐标系统。特征值法被选为主要方法，因为它在处理不同类型地震数据时表现出良好的稳定性和广泛的适用性。

2.1 特征值法在单层介质中的应用

SWSPy中用于测量剪切波分裂的特征值法包括以下步骤（参见图3），针对每个接收器的S波到达，适用于所有地震事件：

  1. 加载数据并进行必要的预处理。
  2. 将数据旋转到LQT（传播、垂直横向、水平横向）坐标系统。
  3. 计算所有可能快方向和延迟时间下的特征值比率 (λ2/λ1)，以获得最佳分裂参数（δt, φ）。
  4. 执行聚类分析，寻找对应最小λ2/λ1的最佳分裂参数。
  5. 如果需要，计算质量度量 QW (Wuestefeld et al., 2010)。
  6. 从剪切波分裂修正的粒子运动中计算S波源的偏振。
  7. 将分裂参数结果从LQT坐标系转换为ZNE坐标系。

2.1.1 预处理

数据预处理是分析中的第一步，主要包括以下内容：

  1. 去趋势和滤波：首先对数据进行去趋势处理，去除数据中的趋势成分，然后进行必要的滤波操作，去除噪声，确保S波信号不受影响。
  2. 上采样和下采样：根据原始数据的采样率和计算效率需求，数据可以进行上采样或下采样。上采样提高了延迟时间（δt）测量的精度，但会增加计算成本；下采样则减少了计算成本，但牺牲了一些精度。
  3. 上采样方法：上采样使用加权平均斜率法，这是一种通过加权平均来提高数据采样率的方法。
  4. 下采样：下采样通过减少数据点数来降低计算量，但会降低δt测量的精度，适用于计算资源有限的情况下。
  5. 仪器响应的去除：如果S波信号超出了仪器的常规响应带宽，需要在预处理阶段移除仪器响应，以确保信号分析的准确性。

2.1.2 旋转到LQT坐标系统

2.1.3 寻找最佳分裂参数

一旦数据旋转完成，可以通过网格搜索来找到最佳的剪切波分裂参数δt和φ，以使数据最线性化（能量在偏振（P）平面内最大化，在空位（A）平面内最小化，见图2）即SC91。

对于每一种可能的δt-φ组合，Q(t)和T(t)首先在QT平面内旋转φ角度，然后分别在时间上前后偏移δt/2。接着，构建Q(t)和T(t)的协方差矩阵，并求解该矩阵的特征值。第一和第二特征值的比率（λ2/λ1）描述了QT平面内粒子运动的线性度，较小的比率表示数据的线性度较高。为了最大化解的稳定性，使用λ2/λ1而不是λ1/λ2。

网格搜索是计算量最大的步骤，计算成本取决于δt和φ的分辨率。为了减少计算成本，SWSPy使用了numba编译器，将执行网格搜索的函数转化为机器代码运行，从而加速计算过程。

2.1.4 多窗口稳定性聚类分析

S波相位的窗口选择对结果的稳定性有显著影响。为了找到最稳定的结果，我们实现了Teanby等（2004）的方法，通过改变窗口的开始和结束时间并进行数据聚类，寻找最稳定的结果。这包括对每个窗口在δt-φ空间中重复进行网格搜索。

多个窗口的示例可以在图5a中看到，用户可以指定窗口的持续时间、开始和结束位置，以及窗口组合的数量（如图4和清单1所示）。对于完全自动化的剪切波分裂分析，必须在处理之前指定这些参数，这与非自动化方法中通常为每个事件单独选择这些参数有所不同。用户通过指定以下参数来控制窗口选择：

这些参数都是相对于S波相位到达时间定义的，这意味着用户必须提供S波相位到达时间的适当近似值。

每个窗口的最佳分裂参数δt和φ通过DBSCAN聚类算法进行聚类（Ester等，1996）。这与Teanby等（2004）的方法不同，因为我们在新的领域中进行聚类，优化处理了φ的周期性特性。聚类领域C的定义如下：

$$C = \left( \delta t \cdot \cos(2\phi), \delta t \cdot \sin(2\phi) \right)$$

其中，δt是归一化的延迟时间，φ是快方向的偏振角度。通过聚类，SWSPy找到每个源接收器对的最优分裂结果，选择具有最小方差的结果。这些方差包括集群内方差和数据方差，其计算方式如下：

$$sigma^2_{\text{cluster},c} = \frac{1}{N_c} \sum_{n=1}^{N_c} \left( (\delta t_n - \bar{\delta t}_c)^2 + (\phi_n - \bar{\phi}_c)^2 \right)$$

$$sigma^2_{\text{data},c} = \left( \sum_{n=1}^{N_c} \frac{1}{\sigma^2_{\delta t,n}} \right)^{-1} + \left( \sum_{n=1}^{N_c} \frac{1}{\sigma^2_{\phi,n}} \right)^{-1}$$

其中， $N_c$ 是集群 $c$ 中的样本数， $\bar{\delta t}_c$ 和 $\bar{\phi}_c$ 分别是集群 $c$ 的平均值。通过这些计算，SWSPy能够准确地聚类分裂参数并获得最稳定的结果。
Image

2.1.5 自动化处理多个地震源和多个接收器

QW值通过以下公式计算：

$$ QW = { \begin{matrix} -(1 - d_{\text{null}}) & \text{当 } d_{\text{null}} < d_{\text{good}} \\ (1 - d_{\text{good}}) & \text{当 } d_{\text{null}} \geq d_{\text{good}} \end{matrix} } $$

其中，dnull和dgood分别表示较差和较好的测量结果。QW值的高低反映了测量结果的质量，理想结果对应QW = 1，而差的结果则接近零。

2.1.6 S波源偏振

一旦获得最佳的剪切波分裂结果，就可以去除剪切波分裂的影响，恢复地震源辐射出的原始S波。初始的S波源偏振可以通过去除各向异性后的S波粒子运动在QT平面中的特征值来获得。S波源偏振是地震分析中一个有用但未被充分利用的参数，因为对于双偶极地震源，它代表了断层滑动的方向。我们在第3.3节中提供了一个关于如何进行源偏振诊断的示例。

2.1.7 从LQT到ZNE坐标系统的旋转

最后，所有的结果，包括最佳的快方向（φ）、各项质量指标以及S波源的偏振，都将从LQT坐标系统转换到ZNE坐标系统（有关所有相关角度的定义见图2）。因此，结果呈现了一个完整的三维结果，这能够为地震分析提供更详细和全面的信息。

2.2 扩展到多层介质的方法

方法概述：

  1. 同时反演两层的方法：
    一些学者（Silver和Savage，1994；Özalaybey和Savage，1994；Wolfe和Silver，1998；Reiss和Rümpker，2017）提出了通过同时反演两层剪切波分裂参数的方法。这些方法通过使用多个地震源信号到达同一个接收器，并结合表观分裂参数与每层实际分裂参数之间的理论关系来反演多层的最佳属性。简单来说，这一过程类似于一维层析成像问题，通过组合多个地震源接收器的数据来反演多层介质的结构。

  2. 各向异性层析成像方法：
    另一种方法是将介质划分成多个箱形域（通常为水平层），每个域都包含完整的各向异性弹性张量，并求解Christoffel方程来计算理论分裂参数。这种方法本质上是一种各向异性层析成像方法，通过反演得到每一层的剪切波分裂参数。与前述同时反演两层的方法相比，这种方法更为稳定，且能够解决更复杂的多层问题，但其一个显著的限制是需要明确指定各向异性层的厚度（Wookey，2012；Hammond等，2014；Kufner等，2023）。

新提出的逐层方法：

2.2.1 逐层方法所需假设

  1. n层将S波分裂n次（Yardley和Crampin，1991；Silver和Savage，1994）。
  2. 每层具有单一的有效各向异性。这意味着该方法只能解析给定层内所有各向异性贡献的整体效果，类似于单层方法。
  3. 最深层（层1）的延迟时间（δt1）必须大于S波的最长主导周期成分（参见Rümpker和Silver（1998b），图1，展示了频率与延迟时间效应的关系）。
  4. 初始表观单层测量中占主导地位的信号来自第一层分裂。这个假设对于大多数情况有效，因为第一层仅在快波和慢波之间分配能量。
  5. 每层的各向异性具有相同的频率依赖性行为，即S波在不同层中的散射行为是相同的。
  6. 每层的快方向（φ1，φ2，...，φn）在QT平面中不平行或正交。如果它们是正交的，则无法区分来自两个层的相位，因为快波和慢波不会进一步分裂，导致其中一个层的结果为零（零结果意味着各向异性不可区分）。

尽管这些条件看起来很严格，但在实际的物理场景中，许多情况可能满足这些假设，因此该方法在实际应用中具有潜力。

2.2.2 两层的分裂方法

多层分裂方法测量每一层的分裂参数（φi，δti），同时使用单层方法测量表观分裂参数（φapp，δtapp），以便量化多层结果相对于单层结果的显著性。这些参数的测量过程如下：

  1. 使用单层方法测量表观分裂参数，对于包含所有S波能量的窗口wininit进行测量（见第2.1节）。
  2. 初始窗口被分成两个子窗口，一个从twininit,start到twininit,start+δtapp，另一个从twininit,start+δtapp到twininit,end，由表观延迟时间δtapp控制。
  3. 对这些子窗口使用特征值法（见第2.1节）分别测量分裂参数，最线性化的结果（最小的λ2/λ1）定义为最浅层（对于两层问题是层2）的最优分裂参数。
  4. 然后对整个窗口wininit中的S波到达时间进行修正，去除层2的分裂影响。
  5. 使用修正后的数据，在整个窗口wininit中再次测量分裂参数，得到最深层（层1）的分裂参数。
  6. 最后，可以确认两层解决方案是否比单层表观解决方案更准确地描述了介质。这里我们定义多层结果更为准确的条件是：（1）更线性（即（λ2/λ1）multi-layer < （λ2/λ1）single-layer）；（2）两层的快方向具有不同的方向，经过不确定性修正后。

多层结果的线性度（λ2/λ1）multi-layer通过对每一层的λ2/λ1值进行求和来量化：

$$ (λ2/λ1)multi-layer = \sum_{n=1}^n \left( \frac{λ2}{λ1} \right)_n, $$

其中，n表示第n层。

2.2.3 扩展到n层

第2.2.2节特地为清晰起见，描述了两层的多层方法。然而，将该方法扩展到n层在理论上是非常简单的。第2到第4步可以通过连续的小窗口重复进行，使用δt2,app, δt3,app, ..., δtn,app来划分每个窗口。然而，实际上可以独立测量的层数是有限的。

2.3 SWSPy使用示例

SWSPy支持从简单的单一源接收器对到多个接收器和多个源的剪切波分裂自动测量。提供了一个简单的例子，展示如何在多个接收器上测量单一源的剪切波分裂，以及如何进行正演模拟，生成具有剪切波分裂的合成信号。

2.3.1 测量地震剪切波分裂

SWSPy使用基于Python类的结构实现（见Listing 1），并广泛利用obspy进行地震数据的输入和输出（Krischer等，2015）。通过传入一个包含所有接收器数据流的obspy数据，创建一个splittingObject，并可以指定窗口和参数搜索空间来执行剪切波分裂分析。

Listing 1展示了如何使用splittingObject来执行剪切波分裂分析：

import swspy, obspy
# 创建剪切波分裂对象：
st = obspy.read(<数据路径>)
splittingObject = swspy.splitting.create_splitting_object(st)
# 指定一些关键参数...
splittingObject.win_S_pick_tolerance = 0.1
splittingObject.overall_win_start_pre_fast_S_pick = 0.3
splittingObject.overall_win_start_post_fast_S_pick = 0.2
splittingObject.max_t_shift_s = 1.0
# 执行分裂分析：
splittingObject.perform_sws_analysis(coord_system=‘‘ZNE’’, sws_method=‘‘EV’’)
# 绘制并保存结果：
# （将splittingObject.sws_result_df保存为csv文件）
splittingObject.plot()
splittingObject.save_result()

2.3.2 正演模拟

SWSPy还支持正演模拟，用于生成通过各向异性介质传播的合成地震波形。Listing 2展示了一个示例，如何生成一个主频为10 Hz的S波合成地震波形，传播通过一个具有60°快方向和δt = 0.5 s的层。

这种正演模拟功能可用于验证SWSPy的性能，并解决反演问题。以下是一个使用正演模拟生成合成地震波形的代码示例：

import swspy
# 创建源时函数：
seismogram_dur_s = 10.0
sampling_rate_hz = 1000.0
st = swspy.splitting.forward_model.create_src_time_func(seismogram_dur_s, sampling_rate_hz)
# 指定层的各向异性参数：
phi_from_N = 60
dt = 0.5
back_azi = 0
event_inclin_angle_at_station = 0
# 应用分裂：
st = swspy.splitting.forward_model.add_splitting(st, phi_from_N, dt, back_azi, event_inclin_angle_at_station)

3 Examples

3.1 简单的冰震示例

在此，我们以冰川上的实际地震为例，展示如何使用SWSPy进行S波分裂分析，特别关注指示可靠测量的关键特征。图5展示了来自南极鲁特福德冰流的单个接收器处的冰震基底粘滑S波到达（Hudson等，2020a；Smith等，2015）。

冰川冰体具有强烈的各向异性结构，结合南极地区的低噪声水平，为S波分裂提供了一个理想的实际案例（Smith等，2017；Harland等，2013；Kufner等，2023）。基底粘滑冰震也提供了理想的例子，因为它们的S波源极化通常受到良好约束，通常与冰流方向大致对齐（相对于北方160°，Smith等，2015），在本例中，通过全波形源机制反演得到了确认（Hudson等，2020a）。
以下是用于量化剪切波分裂结果质量的几个关键特征：

  1. 检查ZNE坐标系中的原始波形与去除分裂后的波形（见图5a）。

    • S波到达波包应位于可能窗口的起始和结束之间（灰色垂直线）。
    • 去除分裂后的波包应比原始数据波包更短。

  2. 最大化和最小化去除分裂后的P和A分量能量（图5b）。

    • P分量与A分量的振幅比表示去除分裂后粒子运动的线性度，较小的λ2/λ1值表示更线性的结果。
    • 对于冰震，λ2/λ1 = 0.033，大部分能量集中在P分量，仅有少量能量到达A分量。

  3. 快慢S波相位到达的时间应不同，去除分裂后时间应对齐（见图5c右面板）。

    • 快慢S波相位应在去除分裂前到达不同的时间，去除分裂后应对齐。

  4. 在东北平面上大致线性的粒子运动（见图5d）。

    • 对于冰震，粒子运动大致线性，源极化大约为从北方165° ± 6°，与冰流方向一致。

  5. 检查聚类分析的稳定性（见图5e）。

    • 聚类样本应有小的不确定性，从而得到稳定的φ和δt解。
    • 如果样本显示明显的非均匀行为，结果可能受到周期跳跃等效应的影响。

  6. 在φ−δt空间内特征值比呈现明显的最小值（见图5f）。

    • 冰震展示了明显的单一全局最小值，最优解由绿色点及其误差条表示。
    • φ−δt空间图有助于检查周期跳跃是否发生。如果周期跳跃占主导，可能会看到多个最小值，并且这些最小值对应的φ值相隔90°。

  7. 测量质量参数λ2/λ1和QW。

    • SWSPy可以输出多个质量参数，λ2/λ1表示结果的线性度。
    • SWSPy还计算Wuestefeld质量因子QW，QW = 1表示良好结果，QW = 0表示较差结果，QW = −1表示零结果。
    • 对于冰震，QW = 0.969，确认了特征值法和交叉相关方法得到的一致结果。
      Image

3.2 远震剪切波分裂

本节展示了SWSPy在远震剪切波分裂中的表现。远震剪切波分裂技术，特别是SKS、PKS和SKKS相位分析，是用于约束上地幔变形模式的常用方法。由于这些相位的近垂直入射角和由P波转化为S波导致的径向极化，使得它们能提供可靠的地幔剪切波分裂测量。

图6展示了2005年2月5日Mw7.1级Celebus海地震的数据，记录于美国加利福尼亚NEE站。使用SHEBA代码进行的剪切波分裂分析发现了SKS和SKKS的剪切波分裂不一致，其中SKS为零结果，而SKKS则表现出明显的剪切波分裂，得到了φ = 74° ± 5°，δt = 1.05 ± 0.07秒的结果。

使用SWSPy重新测量SKKS相位的剪切波分裂，得到了φ = 74.2° ± 14.0°，δt = 1.05 ± 0.175秒，与SHEBA结果一致。我们还能够获取源极化为115° ± 7°，与SHEBA的测量值115°以及观察到的背方位角294°一致，符合SKS径向极化的假设。

这个示例仅展示了一个简单的远震应用，实际中的远震剪切波分裂研究，特别是针对最下层地幔的研究，更加复杂。现代远震剪切波分裂研究采用堆叠和波束形成等预处理方法，有助于在噪声较大的数据集中识别SKS、SKKS和S3KS相位。

3.3 自动化S波分裂分析在火山中多个地震中的应用

Image

多层介质示例

3.4.1 正演模型示例

我们首先在模拟数据上展示新多层分裂方法的表现，然后再将其应用于真实案例。图8展示了一个两层正演模型的结果。剪切波分裂应用于中心频率为10 Hz、源极化方向为0° N的Ricker小波上两次，模拟波传播通过两层介质（φ层1 = 60°，φ层2 = 40°，δt层1 = 0.5秒，δt层2 = 0.2秒）。

图8a-d展示的表观剪切波分裂测量明显无法找到真实结果。然而，φ−δt空间显示表观测量对两个层次都敏感，并且在δt = 0.2秒和δt = 0.5秒处有明显的最小值。第一层表现出较强的分裂信号，因此该层主导了解析。由于周期跳跃的影响，在φ−δt空间中为每一层选择最小值变得具有挑战性。

通过新的层-by-层分裂测量方法，成功解析了两层的各向异性，所有结果都接近真实值，并且大多数结果与真实值在不确定性范围内一致。修正后的波形进一步验证了新方法的性能。

Image
a-d部分：展示了表观的单层测量结果，包括合成的地震波形和分裂参数的估计。
e-h部分：展示了逐层反演的结果，其中g图中的蓝色数据表示仅对第二层进行中间校正后的粒子运动。

3.4.2 冰震示例

成功的多层S波速度各向异性测量在现实世界中的例子较少，这主要是因为测量多层分裂的挑战，而非缺乏多层各向异性介质。然而，冰川冰为多层各向异性提供了一个现实的例子。

图9展示了两层S波分裂结果的水平粒子运动，并与图5中的单层结果进行了对比。结果表明，两层模型比单层模型提供了更好的线性化效果，特征值比λ2/λ1表明，两层结果的线性化程度大约是单层结果的两倍。这证明了两层介质比单层介质更能准确描述观测结果。

修正后的波形进一步验证了我们新方法的表现，结果显示两层介质能够有效地解析冰震中的各向异性，而不是仅仅依靠额外的自由度进行拟合。尽管多层分裂分析提供了更多自由度，但需要小心避免过拟合，尤其是在连续层的快方向一致时。

3.4.3 多层S波分裂的挑战和局限性

因此，尽管多层分裂分析方法展现了潜力，但在实际应用中仍然面临周期跳跃和非唯一解等挑战，使用时需要小心处理。

Image

3.5 SWSPy与其他剪切波分裂软件包的比较

Image

4Disscussion

4.1 优势与局限性

4.2 剪切波分裂的其他应用

Paper link
[Github project address]（https://github.com/TomSHudson/swspy）

转载请注明出处