SDpy:面向稳健地震应力降估计的开源Python软件包

地震应力降连接了观测地震波与断层破裂物理：它既被用于讨论断层强度、破裂效率和震源尺度，也直接影响高频地震动及地震危险性认识。然而，应力降并不是一种可以从波形中直接读取的观测量，而是由拐角频率 fc、地震矩 M0、介质速度和震源模型共同推导的参数。震源谱、传播衰减和场地响应之间的耦合，以及P波、S波和尾波之间的差异，使得同一地震的应力降估计可能因处理方法和模型选择而显著变化。SDpy所面对的核心问题，正是如何把这一长期依赖个人程序、经验参数和分散工作流的分析过程，转化为透明、可配置、可复现并可相互比较的开放软件体系。

    这篇论文的重要性并不在于提出一种新的震源谱理论，而在于将谱拟合法和经验格林函数谱比法统一到同一个Python框架中，并使P波、S波、S波尾波、单窗口与多窗口、单台与多台、单个与多个经验格林函数、Brune与Boatwright等不同方案能够在一致的输入输出体系下运行。由此，研究者不仅可以获得应力降，还可以系统辨别结果究竟受传播衰减、经验格林函数质量、震相选择、频带、窗口还是k系数控制。

    论文对Parkfield重复地震的验证尤其具有方法学价值。P波谱比结果与Abercrombie（2014）的已发表值基本一致；高相关性的经验格林函数与S波尾波多窗口分析显著压缩了拐角频率的离散程度；相反，单谱拟合得到的应力降仅约5–10 MPa，明显低于谱比法的26–70 MPa，直接揭示了简化衰减模型和场地项未充分分离时的系统偏差。与此同时，S波和尾波推导的应力降又普遍高于P波结果，说明 k值及其所代表的震源几何与破裂动力学假设，可能比数值优化本身更强烈地控制最终应力降。

    因此，我推荐将SDpy视为一个“应力降方法学实验平台”，而不是一个能够自动给出唯一真实应力降的黑箱。它的开创性认识是：提高应力降研究的可信度，不能只依靠更复杂的拟合算法，还必须将不同震相、不同震源模型、不同经验格林函数和不同 k值带来的系统差异显式化、可视化并纳入比较。对于区域小震、诱发地震、重复地震和密集台阵研究，这种标准化与可审查性具有直接价值。

应力降通常定义为断层破裂前后静态剪应力之差，它描述断层滑动过程中释放的应力量级，并用于约束断层强度、破裂效率、震源动力学以及地震动特征。由于它连接了观测波形与震源物理，应力降已成为认识地震破裂过程和地震危险性的基本参数。

    但应力降并非直接观测量。基于谱拐角频率的估计必须先从记录中分离震源、传播路径、场地响应和仪器响应。传统单谱拟合通常同时求解震源谱和衰减项，因此拐角频率 fc与品质因子 Q等传播参数之间存在权衡；如果场地效应或频率相关衰减未得到适当描述，震源谱的高频衰减可能被错误归因于较低的拐角频率。

    经验格林函数方法试图绕开这一困难。若一个较小地震与目标地震位置接近、震源机制相似，并由同一台站记录，则两者近似共享传播路径、场地响应和仪器响应。将目标地震频谱除以小地震频谱后，上述非震源项可在一阶近似下消除，从而突出两次地震的相对震源谱。该方法减少了对显式衰减模型的依赖，但要求区域内存在位置足够接近、波形高度相似、震级又明显较小且信噪比充分的经验格林函数事件。

    震相选择同样影响稳定性。P波通常较少受后续震相干扰，波形相对清晰；S波往往具有更高信噪比，但衰减较强；S波尾波由多次散射S波组成，其空间和方位平均效应可降低对台站方位覆盖和目标—EGF微小位置差异的敏感性。因此，论文不是在三者之间预先指定唯一优选震相，而是将其全部纳入统一软件框架。

作者指出，尽管谱拟合、谱比、谱分解、有限断层和地震动反演等多种应力降方法已经存在，实际研究仍缺少一个公开、易用并能够统一不同处理情景的软件包。这一缺口会提高技术门槛，使不同研究组的窗口、频带、震源模型和质量控制难以对照，也不利于结果复现。

SDpy据此采用四层解决思路：

1. 在理论层面同时实现谱拟合法和经验格林函数谱比法；
2. 在数据层面支持P波、S波和S波尾波，以及单窗口和多窗口分析；
3. 在观测配置层面支持单分量或三分量、单台或多台、单个或多个EGF；
4. 在工程层面将数据下载、预处理、多窗谱计算、模型拟合、并行计算、质量控制、绘图和应力降计算封装为模块化工作流。

    这些功能构成论文所称“robust”的主要来源，即通过多方案比较、质量控制和可配置处理降低单一流程的脆弱性，而不是声称任意输入均能自动得到无偏结果。

1 单谱拟合法

    台站记录在时域中可写为震源时间函数 E(t)、传播路径 P(t)和仪器响应 I(t)的卷积；在频域中则为：

    SDpy采用的一般化震源谱为：

其中，Ω0为长周期谱平台，fc为拐角频率，n和 γ控制高频衰减形态。当 n=2,γ=1时为经典Brune模型；当 n=2,γ=2时为Boatwright模型。

    仪器响应去除后，位移谱进一步表示为：

其中 t 为传播时间，论文实现中的Q为频率无关品质因子。SDpy通过最小化观测谱与理论谱在对数域中的平方差确定最优 fc。拐角频率的不确定范围定义为拟合方差不超过最小方差1.05倍时所对应的频率区间。

    应力降由Eshelby圆形裂纹关系计算：

其中 M0为地震矩，β为剪切波速度，k将拐角频率转换为震源半径并依赖所采用的破裂模型。由于 Δσ∝fc^3且 Δσ∝k^−3，拐角频率和 k的小幅变化均会被立方放大。

    论文引用既有研究指出，当真实 fc低于可用频带上限约20%–25%时，高质量数据中的 fc—衰减权衡可明显减弱。这一条件实际上要求频带高端必须显著高于待估拐角频率，因此Nyquist频率、仪器带宽、噪声和高频衰减共同决定可分辨的最大 fc。

2 经验格林函数谱比法

    目标地震与EGF的振幅谱比为：

    若两次地震由同一台站记录且传播路径和仪器响应近似相同，则可简化为：

    将一般化震源谱代入后，谱比由目标与EGF的相对地震矩以及两个拐角频率控制。目标事件的 fcT通过最小化理论谱比和观测谱比在对数域中的差异获得。单个EGF情形下，可将不同台站谱比取均值或中位数后统一拟合；多个台站和多个EGF情形下，则先独立拟合每一条谱比，再以拟合方差的倒数作为权重合并各个 fc测量。

    逆方差加权使拟合较差的测量获得较小权重，并降低最终结果对某个特定台站或EGF筛选阈值的敏感性。不过，这一统计处理不能弥补错误EGF造成的系统偏差；谱比法成立的前提仍是目标事件与EGF具有近似相同的路径、场地和辐射特征。

    作者据Abercrombie（2015）提出三项优选条件：EGF与目标事件空间上非常接近；同台波形具有较高相似性；EGF相对目标事件足够小，同时在拟合频带内仍保持可接受信噪比。在低地震活动区无法满足这些条件时，论文建议谨慎解释谱比结果，或改用对衰减参数施加约束的谱拟合法。

1 数据下载模块

    当本地尚无数据时，SDpy可根据事件元数据从数据中心自动下载SAC格式波形和StationXML仪器响应文件，默认数据中心为IRIS。模块还可根据指定速度模型计算理论P、S到时并写入SAC头段。

2 数据预处理模块

    预处理包括去除仪器响应、计算振幅谱、在缺少 M0 时估算地震矩，以及依据用户设定的信噪比阈值筛选频率点。论文将SNR定义为信号谱平均绝对振幅与噪声谱平均绝对振幅之比，并使用Prieto开发的multitaper程序计算频谱。

    谱窗长度可通过 fixed_window 直接指定；默认情况下，软件依据地震矩按照Abercrombie等（2017）的关系确定窗口长度，若未提供地震矩，则根据地方震级按照Ruhl等（2017）估算。多窗分析中的taper数量和相邻窗口重叠比例分别通过 num_tapers 和 overlap 控制。

3 模型拟合模块

    软件同时实现Nelder–Mead和L-BFGS-B优化算法，并通过 num_workers 指定CPU数量实现并行拟合。相同的拟合框架既适用于单条振幅谱，也适用于共址地震对的谱比。

    源模型可选择Brune模型、Boatwright模型或一般化模型。一般化模型并不同时自由优化 n、γ和 fc，而是在 1.5≤n≤3、1≤γ≤2范围内对固定组合进行系统循环，以减弱谱形参数与 fc的权衡。一次分析中，所选源模型统一应用于整个数据集，以保证事件间比较的一致性；若需研究某一事件的特殊谱形，则需要另行运行。

    用户可采用以Nyquist频率为上限的默认频带，也可指定自定义频率范围。软件只使用同时位于指定频带且满足SNR条件的频率点；通用默认SNR阈值为2，而论文的四个验证实例采用了更严格的SNR >3。

4 可视化和应力降计算模块

    可视化模块将目标事件与EGF波形、时间窗以及谱或谱比拟合结果输出为PDF，并根据单谱拟合、单EGF或多EGF情景自动调整版式。最终计算模块整合预处理、频谱计算、模型拟合和应力降计算，通过导入模块和函数调用即可批量运行。

    所有参数集中存储于JSON格式的 control_file.json。完整流程依次为参数初始化、读取用户设置、计算振幅谱或谱比、估算 fc及其不确定性、代入 M0、k和 β 计算应力降，并将数值结果和图件写入用户指定目录。

图1以流程框图表示三个阶段。左侧“Data Preparation”要求准备目标事件、可用EGF和全部台站的元数据，以及波形、仪器响应和可选速度模型；上方“Parameter Setting”负责指定分析参数；右侧“Stress Drop Estimation”依次完成预处理、P/S/尾波频谱计算、谱或谱比拟合以及应力降计算。

方法1分支先汇总全部台站振幅谱，取中位数或平均谱后进行理论模型拟合。方法2分支计算所有台站和EGF组合的谱比，既可对各谱比单独拟合后按方差加权合并 fc，也可先形成平均或中位谱比再整体拟合。三个分支最终汇入统一的应力降计算模块。图1因此清楚区分了“先堆叠再拟合”和“先逐条拟合再加权”两种统计路径。

作者选取San Andreas Fault Observatory at Depth研究区的五个Parkfield重复地震作为目标事件，使用Northern California Earthquake Data Center高分辨率地震台网的波形。该台网包含13个台站，采样率为250 Hz。五个目标事件震级为 ML2.0–2.2。

    所有P、S到时在实例中均由人工拾取；作者建议在大规模数据处理中使用PhaseNet等自动拾取算法。四个实例分别为：P波单谱拟合、P波多EGF谱比、S波多EGF谱比，以及S波尾波单EGF多窗口谱比。

    四个实例使用的关键模型为：

    所有应力降计算采用与Abercrombie（2014）一致的地震矩，P波取 k=0.32，S波和S波尾波取 k=0.21，剪切波速度取 β=2.3。

图2横坐标为经度，纵坐标为纬度，显示Parkfield地区San Andreas断层、高分辨率地震台网和研究事件。绿色三角形代表HRSN台站，红色五角星代表目标重复地震，黄色五角星代表EGF事件，蓝色方框标出局部放大区。原文图2，第7页。

台站沿San Andreas断层两侧分布。局部插图以200 m比例尺显示目标事件和EGF之间的近距离关系，说明这些事件在空间上具有重复地震和经验格林函数分析所需的共址条件。

图3以事件21402158和EGF11为例，上图为目标事件，下图为EGF。横坐标为相对时间，范围约0–9 s；纵坐标为位移，单位为nm。灰色区域为震前噪声窗Ns，黄色为P波窗，粉色为S波窗，绿色为尾波窗。

P窗和S窗直接截取相应到时附近的波形；尾波窗内部设置五个长度相同且彼此重叠的子窗，不同颜色的水平线表示各子窗位置，相邻窗口重叠50%。目标事件与EGF使用完全对应的窗口，从而保证谱比中的频谱差异主要来自震源项而不是不同的时间截取。

1 实例1：P波单谱拟合

    实例1设置 method=1，由于不使用EGF，必须去除仪器响应，因此 remove_resp=yes。P波谱窗长度为0.5 s，起点位于P到时前0.05 s；所有可用台站谱取中位数，以降低异常台站或异常谱对结果的影响。

    为描述传播衰减，Q被允许在10–500之间变化。五个事件得到的 fc为14.28–15.95 Hz，应力降为5.86–9.13 MPa。该结果显著低于谱比法以及Abercrombie（2014）的结果，作者将其归因于单谱拟合低估 fc，其根源包括用 exp⁡(−πft/Q)描述复杂传播衰减的局限，以及场地效应未被独立消除。

    这一实例并非说明SDpy拟合程序失效，而是显示：当衰减模型不足以代表真实路径和场地效应时，优化算法可以获得数值上良好的拟合，却仍给出偏低的震源拐角频率。

2 实例2：P波多EGF谱比

    作者采用Abercrombie（2014）四个EGF中的三个，另一个因数据无法获取而未使用。每条台站—EGF谱比首先独立拟合 fc及其不确定性，随后将所有测量按逆方差加权；只有拟合方差小于0.005的测量获得非零权重。

    五个目标事件的SDpy结果为 fc=23.28–31.78 Hz、Δσ=26.43–70.44 MPa；Abercrombie（2014）的对应范围为 fc=23.2–30.7 Hz、Δσ=26.0–64.2 MPa。两组结果总体一致，表明SDpy能够复现既有P波EGF分析的主要结果。

3 实例3：S波多EGF谱比

    实例3沿用实例2的窗口长度、EGF集合和逆方差加权方式，仅将P波替换为S波。由于S波高频衰减更强，各台站和EGF得到的 fc离散程度高于P波。五个事件的 fc 为23.08–32.28 Hz，应力降为80.64–267.88 MPa。

    尽管S波 fc 并不总是高于P波，S波应力降却显著更大。主要原因是计算采用 kP=0.32、kS=0.21，而应力降与 k^−3成正比。论文引用动力学震源研究指出，震源几何、破裂方式和破裂速度的差异可使基于固定 k 值得到的应力降相差5–10倍；常见模型中的 k 也可从约0.372变化至0.21。

    因此，P、S结果不一致并不能简单解释为某一震相“错误”。它首先暴露了理论震源模型与真实破裂之间的不匹配，并说明跨震相比较必须同时检验 fcP/fcS和 kP/kS，不能只比较最终的 Δσ。

4 实例4：S波尾波单EGF多窗口谱比

    尾波窗起点设在发震时刻之后三倍S波走时处，适用于本文小于100 km的震中距。考虑到目标事件 ML<3，尾波被划分为五个1 s子窗，相邻子窗重叠50%，形成3 s有效尾波时长。

    作者使用多段匹配滤波方法重新搜索到一个高质量EGF。该事件比目标地震约小一个震级单位，与目标事件的波形互相关系数超过0.95，高于Abercrombie（2014）对同一目标事件所用的EGF。由于这里只使用一个EGF，作者先计算各台站谱比，再取其中位数拟合最终 fc。

    五个事件得到 fc=18.51–22.16 Hz，变化范围明显小于其他实例，应力降为50.50–80.66 MPa。作者认为，尾波的空间平均特性与高波形相似度EGF共同降低了结果离散性；但由于仍采用S波 k=0.21，尾波应力降整体高于P波结果。

图4以事件21402158为例。四个子图横坐标均为频率，单位Hz，并采用对数坐标。图4a纵坐标为位移谱振幅，单位nm/Hz；图4b–d纵坐标为无量纲谱比，同样采用对数坐标。

图4a对应P波单谱拟合。不同细线表示各台站观测谱，粗橙线为堆叠中位谱，绿色虚线为理论源谱，灰线为噪声谱；左下插图显示拟合方差随候选 fc的变化。该事件最优 fc=14.67Hz，不确定度为0.25 Hz。

图4b和4c分别对应P波和S波多EGF谱比。颜色区分台站，线型区分不同EGF，橙色五角星标记每一条谱比独立拟合得到的 fc，插图以柱状方式展示相应方差或拐角频率分布。事件21402158的加权结果分别为P波 24.62±0.25 Hz和S波 23.08±0.32Hz。。

图4d为单EGF尾波谱比。彩色细线表示各台站结果，粗橙线为观测中位谱比，绿色虚线为理论曲线，插图显示方差—拐角频率关系。最优结果为 21.37±0.09 Hz，其形式比P、S多EGF子图更集中，直观支持高质量EGF与尾波分析降低离散性的判断。

对事件21402158，四种方案的应力降分别为5.86、27.69、80.64和64.05 MPa。P波和S波谱比的 fc很接近，但应力降相差近三倍，清楚表明该差异主要来自 k 值而非拐角频率本身。

表2最重要的信息不是哪一列应力降“最大”，而是同一批重复地震在不同方法下可产生一个数量级以内但仍非常显著的系统变化。软件计算的一致性只能保证给定假设下的可复现性，不能消除震源模型、传播模型和 k 值带来的物理非唯一性。

1 谱比法通常优于无约束单谱拟合，但并非无条件成立

    当存在高质量EGF时，谱比法能够同时削弱路径衰减、场地响应和仪器响应的影响，因而比简化衰减模型下的单谱拟合更稳定。本文P波谱比结果与既有研究相符，而单谱拟合明显低估 fc，构成这一判断的直接证据。

    但谱比法的优势依赖于EGF质量。若目标与EGF的空间位置、机制或传播路径差异过大，谱比不能真正消除非震源项；即便拟合曲线光滑，也可能得到有偏结果。作者因此强调EGF选择必须针对具体数据集制定，不存在普遍适用的唯一筛选策略。

2 尾波提供稳定性，但稳定并不等于绝对准确

    多次散射使尾波对辐射方位和小尺度空间错位不那么敏感。本文采用高相关EGF后，尾波 fcf_cfc 的事件间范围缩小到18.51–22.16 Hz，说明其重复性较好。

    然而，尾波应力降仍依赖所选Boatwright谱形、S波 k 值、尾波起始时间、子窗长度和重叠比例。较低的离散性证明测量更稳定，却不能单独证明其绝对值最接近真实应力降。

3 k值是跨震相比较的关键系统误差

    k 将拐角频率转换为震源尺度，其数值隐含断层几何、破裂速度和滑动分布等动力学假设。论文认为P波与S波应力降不一致的主要原因，并非SDpy本身，而是固定 kP/kS与真实破裂过程不匹配。

    作者建议依据观测到的 fcP/fcS比值选择相容的理论 k，而不是机械套用单一模型常数。这一点使论文的意义超出软件发布：它提醒研究者，应力降的跨区域、跨震相和跨目录比较必须首先统一或显式报告震源尺度转换关系。

一、统一双方法框架。 SDpy在同一套输入、质量控制和输出体系中实现谱拟合与谱比，使二者能够针对同一数据直接比较，而不再依赖彼此独立的脚本。

    二、支持多震相和多观测情景。 软件覆盖P波、S波和S波尾波，允许单窗口或多窗口、单台或多台、单分量或三分量、单EGF或多EGF处理。

    三、将多EGF统计合并标准化。 对每条谱比单独拟合并采用逆方差加权，可自动降低低质量测量的贡献，并保留各台站和EGF之间的离散信息。

    四、将尾波多窗分析纳入通用软件。 论文展示了尾波子窗、重叠、台站中位谱比和高相关EGF联合使用的完整实现，并证明该方案可减小重复地震 fc的离散程度。

     五、强调可复现性和可审查性。 JSON控制文件、并行计算、自动输出拟合图及开放源代码，使研究者能够检查窗口、噪声、频带和拟合曲线，而不是只接收最终应力降数值。

首先，论文的实证验证集中于Parkfield一个重复地震序列中的五个 ML2.0–2.2事件。该数据具有高采样率、密集台站和优质EGF条件，因此能够证明软件复现已有分析的能力，但尚不能单独证明在低信噪比、稀疏台网、复杂三维衰减或不同构造环境中均具有相同精度。

    其次，谱拟合所用的频率无关 Q 和指数衰减项是对真实传播过程的简化。实例1本身已经显示，即使允许 Q 在10–500之间变化，源谱和衰减项的权衡仍可能造成系统性 fc低估。

    第三，论文采用“拟合方差不超过最小值5%”定义 fc 不确定性，并在多EGF情况下进行逆方差加权。这是一种操作性明确的拟合不确定度，而不是包含震源模型、 Q结构、场地效应、EGF失配和 k值误差的完整概率后验。

    第四，软件在一次数据集运行中统一采用同一种源模型，有利于事件间一致比较，但也可能掩盖不同事件真实谱形的变化。论文允许通过独立运行考察其他模型或频带，但这不属于默认工作流。

    第五，SDpy的灵活性意味着用户承担更高的方法责任。窗口、频带、SNR阈值、EGF、 Q、 k、 β和源模型均可调整；合理配置能够适应不同数据，不合理配置也可能产生形式上完整但物理上有偏的结果。因此，自动化不能替代波形检查、拟合图审查和多方案敏感性分析。

SDpy建立了一个模块化、可扩展的Python应力降分析框架，将P波、S波和S波尾波的谱拟合与谱比方法，以及不同窗口、台站、分量、EGF和震源谱模型整合到统一工作流中。其数据下载、multitaper频谱计算、并行拟合、逆方差合并、PDF图件输出和JSON参数控制，为应力降研究提供了较完整的可复现基础设施。

    Parkfield重复地震试验表明，P波多EGF谱比能够复现既有研究的拐角频率和应力降；高相关EGF结合S波尾波多窗口分析可降低 fc离散性；单谱拟合则容易受衰减和场地效应影响而低估 fc。P、S和尾波应力降之间的显著差异进一步证明， k 值及其背后的破裂模型是应力降绝对尺度不可忽略的控制因素。

    从地震学发展的角度看，SDpy最重要的贡献不是宣布某一种方法获得了唯一正确的应力降，而是提供了一个可以系统暴露方法依赖性的平台。它使研究者能够将“不同方法为什么得到不同结果”转化为可重复检验的科学问题，并为今后的区域应力降目录、诱发地震源参数研究、重复地震监测以及跨研究结果验证建立共同的软件基础。这是该工作的真正学术价值，也是其值得向地球物理和观测地震学同行推荐的原因。