2026-07-06

地震尾波衰减研究

回顾我那篇论文撰写的历程

姜 - 2022 - 基于单次散射模型的河西走廊地区尾波qc值特征研究

梁 - 2024 - 尾波衰减Qc值在地震预测中的应用研究

主要基于以上两篇硕士学位论文，撰写本文档。

主要使用Visual Studio Code和Antigravity两个编辑器中集成的Claude Sonnet系列的AI模型，以辅助实现本研究中用于数据处理的各类脚本。

1 地震事件

	姜👍	梁
研究区域	河西走廊地区	青藏高原东北缘
时间跨度	2000年6月至 2020年12月	2000年6月至 2023年4月
震级选择	$1.5 \le M \le 5.1$ $( 1.5\le M\le 3$ 和 $3\lt M \le 5.1)$	$M \ge 2.0$
震中距	$\le 100 km$ 内细分三级震中距	$\le 100 km$
震源深度	（未说明）	$1 \le D \le 44 km$
地震记录	14883条地震记录中筛选出5948条地震记录	东部3958次地震西部1252次地震未说明筛选结果

2 台站信息

	姜	梁👍
台站	甘肃省河西走廊区域8个数字固定台网台站、40个主动源流动台、34个喜马拉雅二期流动台	东部：150个喜马拉雅一期流动台站、125个喜马拉雅二期流动台站、59个甘肃地震监测固定台西部：5个甘肃测震台网的固定数字测震台站
频带	（未说明）	宽频地震仪(60s-50Hz、60 s-100Hz和120s-100Hz) 甚宽带地震仪(120s-50Hz)

3 参数设置

	姜👍	梁
波形分量	垂直分量	（未说明）
事件截取	从地震波初至时间前10秒开始，向后截取200秒的地震记录	根据地震发生时刻 $T_0$ 截取地震发生前10s，以及之后300s的地震波形
滤波器	带宽为 $f\pm1/3f$ 的4阶Butterworth带通滤波器	带宽为 $f\pm1/3f$ 的4阶Butterworth带通滤波器
中心频率	1.5Hz、3Hz、6Hz、9Hz、12Hz、18Hz和24Hz	1.5Hz、3Hz、6Hz、9Hz、12Hz、18Hz和24Hz
尾波窗口起始时间	用P波到时及发震时刻确定尾波窗口起始时间 $T_s = 2\sqrt{3}(T_p-T_o)$	对滤波后的数据选取两倍S波走时作为流逝时间窗的开始 $T_s = 2\sqrt{3}(T_p-T_o)$
滑动窗长	动态采样：取滑动窗长为 $\frac{5}{f}$ 秒，步长为滑动窗长的一半 $\frac{5}{2f}$ 秒	动态采样：取滑动窗长为 $\frac{5}{f}$ 秒，步长为滑动窗长的一半 $\frac{5}{2f}$ 秒
流逝时间窗	20s、30s 和40s	20s、30s 和40s
信噪比	$SNR \ge 2$ 和 $SNR \ge 3$ 以进行对比分析	$SNR \ge 3$
最大相关系数	$\le -0.6$ 以及 $\le -0.7$ 以进行对比分析	$\le -0.6$

更为详细的参数设置解释可见朱新运的博士论文《衰减、场地响应等地震波传播相关信息综合研究》

4 流程简介

姜： （数据预处理）
1. 将多种格式的地震波形记录转化成为SAC格式，从中分离出垂直分量地震记录；
2. 挑选出质量较好的波形，通过SAC程序，从地震波初至时间前10秒开始，向后截取200秒的地震记录；
3. 将发震时刻信息等信息写入SAC头段，并手动标出地震波P波到时；
4. 将参考时间由北京时转换为世界时，为Python程序的运行和计算做好准备。
梁：
1. 使用巴特沃夫滤波器对数据进行带通滤波；
2. 滤波后的数据选取两倍S波走时作为流逝时间窗的开始 $T_s$ ，在流逝时间窗上使用滑动窗长进行计算 $Q_c$ ，并计算尾波均方根幅值 $A=(f,t)$ ；
3. 使用最小二乘法拟合与 $t$ 的线性函数包络线斜率 $k=-\pi f/Q_c$ ；
4. 根据不同中心频率和流逝时间窗计算的尾波 $Q_c$ 值，拟合 $Q_c$ 值随频率变化的幂函数关系 $Q_c=Q_0f^{\eta}$ 求取研究区域的 $Q_0$ 值和 $\eta$ 值。
用于计算 $Q_c$ 的diffTW-Qc脚本：
1. SAC数据头段变量要求：
  - evlo: 震源经度 (event longitude)
  - evla: 震源纬度 (event latitude)
  - mag: 震级 (magnitude)
  - o: 发震时刻标记 (origin time marker)
  - b: 记录开始时间 (begin time)
  - t1: P波到时标记 (P-wave arrival time)
  - delta: 采样间隔 (sampling interval)
2. 环境要求：（pyhton环境）
  - ```
    pip install numpy obspy
```

5 处理日志

1记 往日 ~ 2025年10月7日：了解了解地震学？没学什么东西，下载了一点文献，拷贝了点数据。
2记 2025年10月8日 ~ 2025年11月16日：_{此期间基本没什么实际进展，算是玩去了。}
- 基本掌握了地震尾波衰减的波形处理流程；
  - 整理了姜和梁留下来的杂七杂八的数据处理脚本文件，并开始尝试利用AI工具对部分必要的脚本进行重写；
  - 将主动源数据进行入库处理（入库方法介绍已详细撰写在ZDY文件夹下的README.md文档），不然没法提取并使用该数据。入库结果放在ZDY21-24文件夹下，重要的是自动生成的archive.sta归档文件，在后续的提取波形阶段必需该归档文件；
  - 开始进行数据提取，提取结果放置在GetData文件夹下（提取方法介绍已详细撰写在该文件夹下的README.md文档）
    - 该阶段撰写了用于地震事件时间格式转换、地震波形数据提取的脚本；
    - 使用ArcGIS软件从全国地震目录中导出了各个主动源台站100km内的地震事件信息；
    - 通过数据提取脚本，截取得到了40个主动源台站对应100km内的地震波形数据，数据放置在DataSAC300s文件夹下。
  - 撰写了毕业论文《青藏高原东北缘地区震源机制解与应力场分布特征》的开题报告。
3记 2025年11月17日 ~ 2025年11月23日：_{数据预处理阶段，工作文件夹为Preprocessing（预处理介绍已详细撰写在该文件夹下的README.md文档）}。周五下午进行了开题答辩！老师给到的反馈是基本没什么问题了。（2025年12月3日）更正：处于数据提取阶段，工作目录为GatData！
- 撰写用于批量修改SAC数据头段变量的脚本，以得到标记 $o$ 。
4记 2025年11月24日 ~ 2025年11月30日：_{仍处于数据预处理阶段，预计本周的进度会比较缓慢，大部分时间会给到英语六级的复习。}（2025年12月3日）更正：处于数据提取阶段，工作目录为GatData！
- 由于后续计算 $Q_c$ 时还需要输出经纬度和震级，所以还需要在头段变量中写入evlo、evla、mag等信息。目前这些信息都存在筛选过的地震目录中，对于地震目录的格式还需要进一步优化以供这一需求的实现。而对于每一个SAC文件的头段变量的写入方式还在考虑中（是在提取数据时逐个写入，还是后续批量写入？）
  - （2025年11月26日下午） 目前已完成修改用于地震事件时间转换的脚本convert_events_to_DoY.sh。以zdy01-Evnets.txt转换结果文件为例，2001 003 07 37 41 000 40.45 98.49 2.0 7 0分别表示发震年年积日小时分钟秒毫秒纬度经度震级震源深度标识符0且分隔符全部为一个空格；
  - 对于头段变量的修改，我目前的想法是修改getdata_by_station.sh脚本，使用arcfetch提取miniseed数据后不直接批量转换成SAC数据，而尝试一个一个转换成SAC格式并修改头段变量。优点是，在数据提取阶段就已经读取了对应的地震目录信息，此时可以直接写入SAC；
  - 对于这一需求，目前我的计划是使用AI工具来辅助，在提取到miniseed数据后立马转换为SAC数据，并进入sac软件包修改头段变量。（2025年11月27日下午） 才发现mseed2sac软件包（具体操作方式可见该软件包的使用文档）可以直接通过-E选项直接添加当前转换数据的发震年年积日小时分钟秒毫秒纬度经度震源深度，但是对震级信息无法添加，这一信息也许还需考虑如何写入。
- （2025年12月28日）目前采用mseed2sac软件包的-E选项的波形数据转换脚本已完成修改，具体写法见脚本getdata_by_station.sh，后续考虑如何写入震级信息！（还是采取读取文件名匹配的方式来写入？）（2025年12月2日晚上）由于现在引入了自动额外输出偏移后的时间这一特性，写入震级这一需求还是可以通过sac软件包来写入（已实现）。
5记 2025年12月1日 ~ 2025年12月7日：_{处于数据预处理阶段的尾声}，预计本周的进度会比较缓慢，大部分时间会给到英语六级的复习。虽然在此之前说明处于数据处理阶段，工作目录为Preprocessing，但由于地震目录的更换，数据提取GetData不得不进行返工，本周及之前的工作实际处于数据提取阶段！！！成体系地修改了用于该阶段的两个脚本convert_events_to_DoY.sh和getdata_by_station.sh，并更新了该目录下的READEME.md说明文件！！！
- 依据老师的要求，震级范围完全与姜师姐选取的大小一致 $1.5\le M \le 5.1$ ，目前我只有2.0以上的地震目录，这一信息还需和老师信息说明。（2025年12月2日下午） 向老师索取并拿到了研究区内2020年2025年的 $M\le 5.0$ 的地震目录（只需要2021年至2024年全年的地震目录即可），现再次根据各个台站100km内为条件进行地震事件的筛选。为方便后续研究的数据使用（用以求取震源机制解），在完成尾波衰减这一内容后，需要“将 $M \ge 4.0$ 级以上地震波形资料保存好”；
- 当前需要考虑一个问题。在转换世界时如果多减去10s，那么使用-E选项添加o会变得麻烦（需要额外考虑到是否跨60s）。考虑，是否可以在转换世界时的输出中同时包含多减和不减的时间，以供mseed2sac软件包使用？（2025年12月1日晚上） 该需求已得到实现，convert_events_to_DoY.sh脚本在使用-t选项来偏移时间时会自动额外输出偏移后的时间year doy doy_hour doy_min doy_sec jday hour min sec msec evla evlo mag evdp extra。（2025年12月2日下午） 通过该额外输出发震时刻，在使用mseed2sac软件包转换格式时，可以准确当前的文件的发震时刻了；
- 由于地震目录文件内容格式的修改，getdata_by_station.sh脚本当中读取该文件的方式也应注意。（2025年12月2日晚上） 通过sac软件包在该数据提取脚本中集成了在SAC文件头段变量写入震级信息mag的能力；
- （主动源数据问题） 目前发现，zdy13台站编号BA70在2021到2024年并没有任何数据。zdy36台站的miniseed数据转换出来的SAC格式名不是默认的XX.ZDY36.*而是XX.ZDY.*（缺少台站的序号36）这一点需要注意（经查询，姜师姐的主动源数据文件夹36下也存在这一问题）；
- 尝试在sac软件中使用ppk进行P波震相拾取，标记为 $T_1$ ，其余处理见疑问5。同时进行波形筛选，差的波形标记 $T_9$ ，但是目前我对筛选的条件还比较模糊（需和老师进行沟通）。最后针对zdy01台站下部分震相拾取结果进行了 $Q_c$ 计算测试，发现 $Q_c$ 值确实与LTW呈显著的正相关，故提出了疑问7；
6记 2025年12月8日 ~ 2025年12月14日：正式进入数据处理阶段，工作目录位于Preprocessing并进行小范围测试计算 $Q_c$ 值PyQc(byQQL)。主要工作内容为P波震向拾取和使用 DiffTW-Qc 脚本计算 $Q_c$ 值。
- 在进行震相拾取之前，还应进行去均值、去线性趋势和波形尖灭rmean; rtr; taper，而去仪器响应在本研究中是非必要的操作，故不进行transfer处理——使用脚本preprocess_sac.sh即可；
- 在进行震相拾取时，似乎振幅量级在大于 $10^3$ 得情况下比较好，目前建议优选该类量级的波形进行计算Q_c值
  - 在使用脚本计算 $Q_c$ 值时，有个SAC别文件可能会因mseed2sac自动输出的文件名不匹配等问题，造成头段变量中没有写入震级信息，此时需要手动写入；
- （2025年12月11日晚上） 有些事情😊暂时告一段落吧，继续准备英语六级和硕士研究生的学习任务。
7记 2025年12月15日 ~ 2025年12月21日：继续数据处理阶段（P波震相拾取），工作目录位于Preprocessing。
- 在利用ppk拾取P波到时时，请注意，截取的波形总时长和P波到时的时间倍数关系！比如截取200s的波形，如果P波到时在90s，则无法计算尾波，因为计算 $2\sqrt{3}(T_p-T_o)$ 的时长已经超出；
- 开始思考如何拟合 $Q_c$ 对频率的依赖关系式，详细内容可见后续的 **如何拟合频率依赖关系 一节；
- 本周在进行P波震向拾取时，从zdy12的数据开始，将会采用更加激进的拾取策略（只关注该次地震波形的总体质量）；
- （2025年12月18日下午） 截至目前，只完成了14个主动源台站数据的震相拾取，但是考虑到后续还有数千个波形数据仍需拾取，可以在后续工作中尝试使用机器学习、深度学习等方式方法来批量自动进行震向拾取；
- （2025年12月19日下午） 发现zdy16这一台站下截取的数据，绝大部分的波形质量都比较差，怯于拾取震相并使用。而与之形成鲜明对比的是zdy17的数据，以至于几乎每一个波形数据的质量都非常好😯；
- （2025年12月21日晚上）：截至目前，已经处理了21个主动源台站的数据。
8记 2025年12月22日 ~ 2025年12月28日：继续数据处理阶段（P波震相拾取），工作目录位于Preprocessing。
- （2025年12月23日下午） 主动源zdy35到zdy40这6个台站的数据过多（基本上每个台站被截取出来的记录都有1400条左右），预计这一周应该只能标记完这些内容；
- （2025年12月24日下午） 对于一些截取出来的200s波形数据，在同一条波形中有两个甚至以上个数的地震记录，比如在zdy35的2022年第7和第8日就有大量的多地震波记录，这一情况下应该如何处理？同疑问8；
- （2025年12月25日下午） 截至目前，对于主动源数据的P波震相拾取的工作内容已完成37/40，后续3个台站还有大约3400条波形记录资料，只能预计在本周结束之前完成处理；
- （2025年12月27日下午） 截至目前，2021年至2024年40个主动源台站的 $M \ge 1.5$ 地震波形资料已全部完成P波拾取，需要注意的是zdy13这一台站在此期间没有任何的数据资料！
9记 2025年12月19日 ~ 2026年1月4日：计算地震尾波的 $Q_c$ 值，工作目录位于PyQc(byQQL)。
- 本周预计会和老师讨论一下接下来数据处理以及测试数据的 $Q_c$ 值计算结果！后续工作预计会接续姜师姐已完成的部分；
- （2025年12月30日上午） 由于暂时还未找到姜师姐之前计算得到的 $Q_c$ 结果文件，故需要根据她留下来的波形数据重新计算 $Q_c$ 值。与之前我处理的波形数据不同的是，师姐没有对SAC数据写入elvo、evla、mag头段变量，并且其对P波到时采用t0标记。在完成对此类不同情况的PyQc代码修改后，已成功计算出遗留下来的地震波形数据 $Q_c$ 值，筛选条件依然为 $co \le -0.6$ 并且 $SNR \ge 3$ ；
- （2025年12月31日下午） 准备跨年，2026，新年快乐！
10记 2026年1月4日 ~ 2026年1月11日：_{开始使用最小二乘法进行拟合。} 由于在此期间发现了诸多用于计算 $Q_c$ 值脚本的小问题，所以需要返工计算地震尾波的 $Q_c$ 值，工作目录位于PyQc(byQQL)以及PyQc(byJXX)。
- （2026年1月4日上午） 经多次计算测试，发现本次P波震相拾取的 $Qc$ 值结果偏低过多，需要审视一下往日的拾取策略，并对之前拾取的波形数据进行审查！目前记录一下可能导致该问题的：
  - 之前7记中的更加激进的拾取策略是不恰当的，必须理应仅保留P波到时明显的波形；
  - 针对200s内有多个地震波形记录的情况，理应舍弃该次波形记录；
- *（2026年1月4日下午）*经过仔细的代码审查发现，我复写的diffTW-Qc.py脚本（依据往日的脚本进行重构复现）中，关于计算SNR的部分算法不适用于当前主动源台站裁切出来的波形数据！虽然我之前对该脚本进行了测试，结果显示该脚本是可信的，但是用于测试的数据是往日数据的时间格式，没有发现该问题。现如今我自己处理数据的时间格式发生了变化，所以很大可能是由于该问题导致了我上午发现的 $Q_c$ 值 偏低过多的问题：
  - 数据问题：往日数据裁切后的发震时刻o为0秒，则数据开始时间点b则为-10秒，而如今我裁剪出来的数据发震时刻o默认为10秒，则数据开始b为0秒；
  - 算法问题：往日在其数据时间格式下，会减去一个开始时间b的-10（即加上10秒），但如今我的开始时间b是0秒，无法得到尾波开始时间点；
  - 解决方案：目前尝试对计算 $SNR$ 和 $Q_c$ 的算法进行修改，但目前还无法保证会对后续计算产生未知的其他影响（理论上不会）。
    1. 撰写新的脚本对SAC数据进行时间上的转换，使之与姜师姐处理得到的数据时间格式一致；
    2. 对当前用于计算 $SNR$ 和 $Q_c$ 的算法部分进行修改，使之能够匹配上尾波的开始时间节点。
```
### 请注意，此处展示的对比仅是往日代码的修改示例
########仅需要关注式中o和b的加减关系#########
# 原有用于计算SNR的部分代码（相对于o）
te1 = 2 * (tr11.stats.sac.t0 - tr11.stats.sac.o) - tr11.stats.sac.b + 20 - 5 # 疑问9：此处标记到时使用错误
# 修改时间基准后的计算SNR代码（相对于b）
te1 = (tr11.stats.sac.o - tr11.stats.sac.b) + 2 * (tr11.stats.sac.t0 - tr11.stats.sac.o) + 20 - 5 # 疑问9：此处标记到时使用错误
# 应用到当前，现在o的值为10而b是0，所以得修改相对时间
coda_end_time = (2 * S_WAVE_FACTOR * (trace_for_analysis.stats.sac.t1 - trace_for_analysis.stats.sac.o) + trace_for_analysis.stats.sac.o + CODA_WINDOW_LENGTH - CODA_TAIL_OFFSET) # 疑问9：正确计算方法

# 原有用于计算Qc的部分代码（相对于o）
ts2_i = 2 * 1.732 * (t0 - o) - b + Lslidewindow * (i - 1)
# 修改时间基准后的计算Qc代码（相对于b）
ts2_i = (o - b) + 2 * 1.732 * (t0 - o) + Lslidewindow * (i - 1)
# 应用到当前，现在o的值为10而b是0，所以得修改相对时间
window_start_time = (2 * S_WAVE_FACTOR *(trace_for_analysis.stats.sac.t1 - trace_for_analysis.stats.sac.o) + trace_for_analysis.stats.sac.o + slide_window_length * window_index)
```
      请注意，此处有两个与往日代码算法不一致的地方：
      - 在计算 $SNR$ 时，往日直接采用P波到时的两倍，而在当前则是2倍S波到时（2倍根号3）；
      - 在计算 $Q_c$ 时，往日代码采用i-1进行索引，而当前则直接使用全部索引（不减1）。
      后续再有问题的情况下，需要回溯这一部分的计算实现，此处标记为疑问9。
- （2026年1月5日上午） 通过使用上述示例修改得到的diff-Qc-v3.1.py脚本，重新计算了我处理的40个主动源台站数据的 $Q_c$ 值；
- （2026年1月6日下午） 需要修正往日代码中对于计算信噪比 $SNR$ 时，最后5秒尾波的python代码的走时错误！问题说明详情请查看疑问9。
- （2026年1月9日上午） 为确保后续研究工作中对 $Q_c$ 值计算方法的完全理解，于今日起开始对用于计算脚本的新一轮审查。对于重构的计算脚本diffTW-Qc-old.py，利用往日的数据复现数据结果，显示完全一致，通过初步审查。以下是需要特别注意的点：
  - 师兄师姐在计算震前噪声前5秒数据段的时候，理应是使用P波到时前5秒，但是在实际代码计算中却使用的是发震时刻o的前5秒。此处需要修改成使用P波到时前5秒，但实际均值计算中的影响不会特别大，可能误差范围在可接受范围内，但 $Q_c$ 均值会比之前的结果更小（ $SNR$ 小得多）；✔已完成修正到diffTW-Qc.py！
  - 师兄师姐在计算 $SNR$ 信噪比最后5秒的python代码处应当采用当前则是2倍S波到时（2倍根号3）这一原则，师兄师姐他们在此处有一个错误，需要注意！正确地采用**2倍根号3的P波走时(即2倍S波走时)**来计算信噪比，会比错误地采用~2倍P波走时~得到的符合筛选条件的记录少。这是因为错误地采用~2倍P波走时~时，用于计算信噪比的波记录大概率并不在尾波时间窗之内，导致这样得到的 $SNR$ 偏大，即同等筛选条件下（ $SNR\ge 3$ ）的记录会偏多；✔已完成修正到diffTW-Qc.py！
  - 师兄师姐往日的数据时间参考格式与我现在处理得到的数据参考格式不同，计算脚本需要小范围修改时间参考部分的代码。✔已完成修正到diffTW-Qc-v3.1.py！
- （2026年1月9日晚上） 目前已完成对所有数据的 $Q_c$ 值计算，修改前后的 $Q_c$ 均值变化不大，但这几天的工作非常重要，修正了往日的部分不当的地方！
- 预计本周结束前完成对所有数据的 $Q_c$ 值计算，并在下周开始进行拟合。
11记 2026年1月12日 ~ 2026年1月18日：上周调整并完成了对所有数据的 $Q_c$ 值计算（存储目录位于PyQc/PyQc(byQQL)以及PyQc/PyQc(byJXX)），预计继续在该工作目录下完成对 $Q_c$ 数据的的拟合。
- （2026年1月12日上午） 通过整理所有数据的当前 $Q_c$ 计算值来看，我处理并得到的40个主动源数据的 $Q_c$ 值在3个时间窗下的 $Q_c$ 均值有以下情况：
  - 低频值稍微偏低（比如中心频率为1.5Hz时，通过师姐的数据计算得到三个时窗下的 $Q_c$ 均值分别约为95、118、147，而我的数据约为89、122、152）；
  - 中间频率的 $Q_c$ 均值基本一致；
  - 高频值基本一致，但与前人研究（参考姜和梁）的结果（高频下的 $Q_c$ 均值在三个时窗下都在2100左右）相比，而我计算得到均值仅约1800左右；
  - 小记： 虽存在多处数据的情况不一致，但在考虑前人在计算信噪比 $SNR$ 时出现的错误，仍应坚持自己的数据计算结果。
- （2026年1月12日下午） 准备拟合各流逝时间窗下尾波 $Q_c$ 值随频率变化的幂函数关系：
  - 拟合前的数据是否还需要清洗，比如除去三倍标准差的 $Q_c$ 值？
  - 利用脚本Qc-fitting.py已完成Q_c值关于频率依赖的函数拟合，且结果比较合理的；
    - 20s： $Q_c=\left(75.9\pm11.5\right)f^{\left(1.04\pm0.06\right)}$ ， $R^2=0.982$ ；
    - 30s： $Q_c=\left(105.0\pm16.9\right)f^{\left(0.93\pm0.07\right)}$ ， $R^2=0.976$ ；
    - 40s： $Q_c=\left(132.0\pm20.0\right)f^{\left(0.87\pm0.06\right)}$ ， $R^2=0.974$ ；
- *（2026年1月14日上午）*学习了如何通过地震射线路径使用克里金插值法来获取研究区二维 $Q_c$ 值空间分布图，详见 *** 如何使用克里金法进行空间插值 一节。目前我对改内容的实施步骤比较清楚，但认为其比较繁琐，特别是关于得到数千条地震记录的地震射线，并由此线转点要素。
  - 对于地震射线的获取，首先需要台站和震中的经纬度信息。目前我处理的数据关于此信息在结果和过程SAC文件中均有写入，而姜处理得到的SAC数据中暂无震中的经纬度信息，后续可能还需了解梁是如何实现这一部分的内容；
  - 对于克里金插值法各个参数的设置，梁虽有提到，但无对其的详细解释。
- （2026年1月15日下午） 工作目录为RayPaths，针对我处理得到的数据撰写用于获取台站-震源射线路径的经纬度脚本文件。以 $Q_c$ 结果文件QcResult.txt（包含震源信息）为主表，循环每一行，从文件名中解析出“台站代码”，然后去StationInfo.csv（台站信息）中查找对应的经纬度，最后将“台站坐标”与“震源坐标”拼接到同一行。
- （2026年1月16日上午） 对梁的处理过的数据进行查看，尝试学习其是如何使用克里金插值法来获取二维 $Q_c$ 值空间分布图的。很遗憾，并未通过该方式得到比较详细的操作步骤，仅找到了其使用的地震射线数据，而其论文中的插值结果栅格图没有找到。
12记 2026年1月19日 ~ 2026年1月25日：继续探索 $Q_c$ 值二维图像化的操作办法。
- （2026年1月19日下午） 为了解决直接使用点插值出来的结果有射线路径的迹象，应该创建渔网来对同一区域内的采样点进行均值化，从而消除了射线路径点的大权重反演影响。详细操作办法请见 *** 如何使用克里金插值法进行空间插值 一节，当前时间节点下，我对二维空间化 $Q_c$ 值的操作步骤已经基本掌握，但是似乎还有相关数据量的问题。
- （2026年1月20日下午） 根据目前已有的数据情况来看，姜处理后的SAC数据中没有台站和震中的经纬度信息，我也没有其使用的地震事件目录，这对后续二维图像化Q_c值工作是一个问题。当前我仅对2021至2024年主动源40个台站的数据有处理结果，显然这是不足够的。为了弥补当前数据量的不足，开始对固定台站的数据进行处理：
  - 时间范围为2021年1月1日至2024年12月31日；
  - 研究区内的台站数量大约有100个，但在我拿到的地震目录以及筛选100km内的地震事件两个限制下，有地震数据的台站为59个（目前已完成对各个台站的地震事件筛选）；
  - 目前我对该数据处理的流程是比较熟悉的，唯一有所顾虑的数据存储问题。当前存储条件下，固定台4年的总数据（SEED格式）大小为4TB左右，但从该数据格式转换到SAC数据格式后，需要更多的存储空间；
- （2026年1月21日下午） 开始准备从固定台站进行对应时段波形数据的提取，工作目录位于GDT/：
  - 由于该类台站数据的提取方式有所不同，故新建了用于该SEED数据提取的地震目录格式转换脚本convert_format.sh，转换输出北京时间、GMT时间、用于裁剪记录的秒数、震级和震源深度；
  - 固定台站的原始数据位于/2021、/2022、/2023、/2024四个文件夹下，数据格式为SEED格式。_{预计使用mseed2sac工具将其转换为SAC格式，此时不使用rdseed工具原因为“不需要仪器响应文件，且我对mseed2sac软件包更加熟悉”} （2026年1月21日下午） 决定使用rdseed工具包，为之前的偷懒行为感到尴尬😜；
  - 针对该SEED数据使用时mseed2sac软件包的操作建议如下：
    - 亦如之前处理主动源数据，在转换数据时将地震事件进行写入头段变量（使用-E选项），请注意在此处使用的是GMT格林威治标准时间的年年积日时分秒/经纬度；
    - 使用mseed2sac工具成功转换成SAC数据后，只保留BNZ垂直通道，删除其他两个通道的SAC数据（此处需要考虑的时如何识别处转换后的SAC文件，因为默认命名与匹配的时间有较大的差距），以减少数据存储的压力；
    - 立即使用sac软件包对该SAC数据进行必要的头段变量的修改，包括但不仅限于发震时刻 $o$ （GMT格林威治标准时间）、震级 $mag$ 的修改；
    - 然后立即使用sac软件包对数据进行裁剪，使用cut o -10 200命令读取该SAC文件，使用ch allt (0 - &1,b&) iztype IB修改数据的相对时间，使用rmean; rtr; taper命令对数据进行简单的预处理，使用w "$newname"对数据进行重命名写入。
  - 在AI工具的协助下，撰写了用于从SEED数据提取并处理的 _{python脚本convert_seed_to_sac.py} 和 _{bash脚本convert_seed_to_sac.sh} （已被移除或修改）：
    - 完全实现了上条记录中的操作建议，比如使用-E选项、仅保留BHZ垂直通道、头段变量的修改以及数据的裁剪和简单预处理；
    - 由于SEED数据内部的数据可能会出现各类问题：
      - 问题：没有任何数据，或者有成百上千个甚至数万个多个小段的数据；
      - _{说明：目前针对多个小段的数据（非完整的数据），我的想法是直接舍去，此类数据处理起来似乎比较麻烦；}
      - _{解决方案：python脚本则会直接跳过该类数据（设置了超时阈值），而bash脚本会一直卡在这一进程（不建议使用，除非保证数据完整）。}
    - 经过测试，还是决定使用rdseed工具来对数据进行转换，缘由是该工具直接支持对SEED格式下内部多段数据的合并，并且同等情况下比mseed2sac工具的读写速度更快，开始对先前的脚本进行修改bash脚本extract_seed_to_sac.sh，预计移除python脚本：
      - rdseed软件包的使用为rdseed -f data.seed -d -o 1 -b 5000000000，使用-b参数将缓冲区设大的情况下，rdseed直接合成完整的SAC数据；
      - 如果多段数据之间有实际的时间中断（例如机器重启、数采丢包），rdseed不会自动补零或插值来强行合并它们，则需要使用sac软件包先读取所有该通道的SAC数据，再合并merge gap zero overlap average，最后写入到SAC文件（该情况下，输出的SAC数据文件名在毫秒处为0000，以此与正常数据有一定的区分度）；
      - 最后对唯一的垂直分量BHZ数据（不论是rdseed软件包默认合成的SAC，还是我使用sac软件包手动均值融合的SAC文件）进行写入头段变量（此处需包括地震事件的发震GMT时间 $o$ 、经纬度 $evla/evlo$ 和震级 $mag$ ）和数据基本处理（裁剪-10 ~ 200的波形、时区校正ChnHdr ALLT v和去均值、去线性趋势以及波形尖灭rmean; rtr; taper）。
- （2026年1月22日上午） 当前运行用于提取固定台数据的脚本extract_seed_to_sac.sh大概需要5 ~ 6小时，所幸经过一晚上的运行，GDT/DataSAC文件下的数据已全部提取（预计使用59个台站，约20000条地震记录），文件大小约为1.8GB。从今日起再次开启对固定台数据的P波初动标记：
  - GS.BYT台站附近是很多地震的震源中心吗？多数波形记录在发震后几秒就观测到了P波初至；
  - 由于筛选的是100km范围的地震记录，仅标记保留 P波初动在35秒之前（一般在30秒之前） 的波形数据；
  - （2026年1月23日下午） 截至目前，还剩下约5400条地震波形记录未进行标记，预计在本周结束之前完成该项任务，并使用diffTW-Qc-v3.1.py脚本计算固定台站所有数据的 $Q_c$ 值。
- （2026年1月25日下午） 本周周日下午4点左右，完成对固定台站数据的P波标记，最终选取标记了7000多条地震记录用于计算Q_c值，总结一下目前的研究情况：
  - 截至目前，我拥有三类数据用于青藏高原东北缘地区的地震尾波衰减研究，分别是我处理的40个主动源台站数据（2021 ~ 2024）、59个固定台站数据（2021 ~ 2024）和姜处理得到很多个台站的数据（2000 ~ 2020，实际我仅拿到了约4000多条可用于计算的数据）；
  - 对于尾波衰减的研究，依据上述三类数据，按计划实现两种结果：
    1. 对该研究区计算出不同时间窗下各个频率的 $Q_c$ 平均值，并使用最小二乘法拟合 $Q_c$ 值对频率的依赖关系式；
    2. 根据上万条符合要求的地震记录计算得到的 $Q_c$ 值，按地震射线路径，使用克里金插值法得到研究区的 $Q_c$ 值空间分布图。
  - 针对上述两种研究结果，预计在下一周完成！
13记 （2026年1月26日 ~ 2026年2月1日） 预计这一周为本学期的倒数第二周，并正式进入研究结果阶段（2021至2024年的地震波形数据已全部处理完毕）：
1. 已全部计算✔不同时间窗下各个频率的 $Q_c$ 平均值
2. 已全部拟合✔ $Q_c$ 值对频率的依赖关系式
3. 正在探索……📊研究区的 $Q_c$ 值空间分布图
- （2026年1月28日上午） 本周前两日开始使用克里金插值法来获取研究区的 $Q_c$ 值空间分布图，但我觉得效果较差，插值出来的空间分布图的地震射线路径还是比较明显（即使已经尝试通过渔网来均值化局部区域的 $Q_c$ 值）。根据老师的建议，将研究区缩小到主动源台站覆盖的范围。
- （2026年1月29日下午） 在AI工具的帮助下，“如何使用普通克里金进行空间插值”的参数设置有了进一步的理解，相关内容预计将写入***** 如何使用克里金法进行空间插值**一节的最后。
- （2026年1月30日下午） 为方便后续克里金插值时使用三个流逝时窗下不同频率的Q_c值，我通过使用ArcGIS的模型构建器创建了一个直接生成用于插值的一键化工具，主要简化了一下工作内容：
  - 按频率筛选数据（每个时窗下都有7个频率需要进行处理）；
  - 生成射线路径（XY转线并投影到WGS84-UTM47N）；
  - 空间采样与降噪（沿线生成点，然后使用特定大小的渔网平均，最后要素转点）。
14记 （2026年2月2日 ~ 2026年2月8日）本周为本学期最后一周，预计于2月7日回家，因而小论文发表任务只能被推迟到春节之后。BGM: Never Letting Go
- （2026年2月6日） 明天就出发回家，今天将抽出时间对目前已完成的工作内容进行整理。
  1. 所有台站的尾波 $Q_c$ 值结果文件均只保存在D:\Desktop\matterMaster\Qiulin,Qin\文件夹当中，其中4-地震尾波衰减文件夹存放了关于 $Q_c$ 值拟合结果和所有台站的地震射线路径及其对应脚本；
  2. 预计所有的波形数据均只存放在桌面硬盘，由于该硬盘不方便携带，预计不会被我带回；
  3. 应用于GIS程序的工程文件夹均存放在D:\Desktop\matterMaster\Qiulin,Qin\当中的研究区工程文件夹当中。目前该文件夹存有插值、地震目录选择等三个工程；

15记 （2026年3月6日 ~ 2026年3月8日） 本周为研二秋季学期的第零周，春节一个月约等于什么也没干。
16记 （2026年3月9日 ~ 2026年3月15日） 本周为研二秋季学期的第一周
17记 （2026年3月16日 ~ 2026年3月22日） 本周为研二秋季学期的第二周，希望这周我的小论文初稿能有较大的进展。BGM: Taki - Yuma.Play
- 通过AI（主要是Gemini）检索了解到，关于地震尾波成像这块最好还是使用 $Q_c^{-1}$ 这一可线性叠加的值来进行研究比较合理，详细见疑问10和11；
- 对部分以前的工作进行了修改，例如如何拟合频率依赖关系这一部分的工作内容中，旧有的使用 $Q_c$ 这一错误的数据计算部分重写成了 $Q_c^{-1}$ 值的倒数利用，包括使用其计算得到的均值和标准差；
- 准备修改用于流程化计算 $Q_c^{-1}$ 的ArcGIS建模工具箱；
18记 （2026年3月23日 ~ 2026年3月29日） 本周为研二秋季学期的第三周，预计本周完成小论文的初稿撰写。
- （2026年3月23日下午） 梳理一下目前关于地震尾波衰减研究中 $Q_c$ 和 $Q_c^{-1}$ 的使用进行简单说明：
  1. 对于20-s LTW、30-s LTW和40-s LTW下各中心频率点所对应的尾波 $Q_c$ 值，目前该均值指的是对数均值，后续建议加上样本量的统计；
    参考文献 1：直接支持“对数平均”在尾波 $Q$ 参数中的应用
    - 文献信息： Zhang, F., & Papageorgiou, A. S. (2010). Attenuation characteristics of Taiwan: Estimation of coda Q, S-wave Q, scattering Q, intrinsic Q, and scattering coefficient. Seismological Research Letters, 81(5), 741-751.
    - 支持理由： 该文章在计算台湾地区平均的尾波 $Q_0$ 和频率依赖指数 $n$ （即 $η$ ）时，明确指出他们对不同台站的参数使用了对数平均，以获得区域的代表值。
    - 文献原句摘录：“Considering the logarithmic average of the estimated values of the parameters $Q_0$ and $n$ of the 14 stations listed in Table 1, we obtain $Q_0=134.41$ and $n=0.72$ .” (参考该文的结论与对比部分)
    参考文献 2：支持尾波特征计算中“对数平均等效于几何平均”及抑制极值的作用
    - 文献信息： Grendas, I., et al. (2022). Can site effects be estimated with respect to a distant reference station? Performance of the spectral factorization of coda waves. Geophysical Journal International, 230(1), 585-613.
    - 支持理由： 该文章在处理尾波能量和场地放大系数时，详细阐述了为什么使用对数平均（几何平均），并在文中直接将二者画等号。这是因为衰减/放大效应是乘性效应（Multiplicative effect），乘性效应的统计特征服从对数正态分布，因此必须用几何/对数平均。
    - 文献原句摘录：“The geometric mean value of the noise energy… is considered as the average noise energy level…” 并在图表说明中明确标注：“…from the corresponding logarithmic average (geometric mean).”
  2. 对于 $Q_c$ 值的频率依赖公式，目前拟合算法中的均值也同样采用的是对数均值；（参考文献同上）
  3. 对于 $Q_c$ 值的空间分布图像，目前在ArcGIS中的自建工具中使用的是 $Q_c^{-1}$ 的普通克里金插值并计算倒数后的 $Q_c$ 分布图；
  4. 对于 $Q_c$ 值的时间变化曲线图，目前 $Q_c$ 的算数平均数和对数平均数都有绘制，后续再考虑使用哪一种。
- （2026年3月25日） 从今日起，开始进行论文的绘图工作（\Qiulin,Qin\研究区工程\制图工程）。目前小论文中的所有图还没有进度，关于使用ArcGIS、Python还是GMT有待考虑：
  1. 波形处理与参数拟合示例图：通过Gemini 3.1 pro撰写绘图代码实现计划书，使用VS Code的 Claude haiku 4.5依据计划书撰写实现代码。目前实现的效果我很满意，后续可在此基础上修改（使用Adobe Illustrator进行微调或排版）；
  2. 频率依赖性拟合图：目前使用Python进行绘制（实现代码有两版，可供后续选择或修改）；
  3. 研究区概况以及台站分布图：目前考虑使用ArcGIS出图，但是还未获取完相关绘图数据（已获取高程DEM栅格数据、全国省市县、站点位置、全国断层矢量数据）；
  4. 时间序列图：目前通过Excel对年际对数均值进行了简单的统计和初步绘制，不是特别美观，后续考虑通过Python再次绘制，风格尽量与前面的波形图和拟合图保持一致。
19记 （2026年3月30日 ~ 2026年4月5日） 无
20记 （2026年4月6日 ~ 2026年4月12日） 本周为研二秋季学期的第五周，说句实话，这学期过去这么四周，这份日志我断断续续的已经不知道如何写下去了。
- （2026年4月8日） 今日在Google Gemini的协助下，对研究区内的地质年代资料进行了初步的分析。实在是不知道这份小论文还得多少时间才能基本完成，虽然我早想着做完了。😊
21记 （2026年4月13日 ~ 2026年5月10日） 这一个月的为本学期的第6、7、8、9周，啥也没干，正在尝试将初稿的“结果”和“讨论”部分分开论述。到此，这份日志我是真有些写不下去。
22记 （2026年5月11日 ~ 2026年5月17日） 本周为第10周，可以说是期中了。为了追赶学习进度，不得已将本周末设为Deadline，希望我能按时完成。
23记 （2026年5月18日 ~ 2026年5月24日） 本周为第11周，上周的论文撰写任务基本完成😴，目前新增修改并尝试投稿的Deadline——下周结束之前（本月底）。
24记 （2026年5月25日 ~ 2026年5月31日） 本周为第12周，虽然这几天写论文的效率比较高，但是我觉得越来越不喜欢干这个事情了。
- （2026年5月27日） 目前小论文基本成型，还需要改进的点：
  1. 正文部分所参考的文献过于久远（例如，关于河西走廊的地质地形地震活动的介绍），建议更新为近几年的论文；
  2. 目前预计投稿中文期刊，先将图中标注改为中文；
  3. 关于研究区概况图中的断层不能标注为“活断层”，统一图中标注的地震震级；
25记 （2026年6月1日 ~ 2026年6月30日） 原谅我自己中间跳过了好几周的记录，不过期间我也是在老师的同意下投稿到中文期刊了，希望能一次录用吧，我不想在学习了，我心好累。
- （2026年6月30日） 我真TM傻逼，把Windows系统下我的中间处理数据和结果绘图脚本以及图件给删除并清空了。
  1. 当天下午想通过DiskGenius程序恢复文件，最终我还是天真了，败在了固态硬盘的TRIM数据清理指令协议。当我点击清空回收站时，我才意识到大事不妙，立即停止清理，但我最重要的Qiulin,Qin文件夹还是被干掉了，文件的编码块全部变成了00 00 00 00。
  2. 比较幸运的是，之前所有的脚本都是通过VS Code编写并修改的，这个程序自动保留了备份，只不过需要我手动召回并辨认脚本的功能（文件名并非我设置的，而是程序自动编码的名称）。如果没有这个功能，我之前的工作全部化为灰烬，我不知道该怎么办了。
  3. 除了脚本的手动恢复，后续这几个月我还得重新制作中间的数据，包括但不仅限于 $Q_c$ 的ArcGIS插值数据、统计数据以及地震事件和台站数据。
  4. 经历这件事情之后，我已经开始害怕将未来的工作数据仅保留一份在某一个存储空间上了。
  5. 再次说明一下主动源台站的数据问题（台站zdy13、zdy26、zdy36），zdy13 台站编号BA70在2021到2024年并没有任何数据，zdy26 台站的miniseed数据转换出来的SAC格式名不是默认的XX.ZDY26.*而是XX.BA49.*，zdy36 台站的miniseed数据转换出来的SAC格式名不是默认的XX.ZDY36.*而是XX.ZDY.*。
26记 （2026年7月1日 ~ 2026年7月5日） 经过几天中午不睡午觉的努力，画图的绝大部分代码都通过VS Code的History文件夹搜索并找回了，但是中间用于插值 $Q_c$ 的数据文件全部都无法恢复了。只能说我得给VS Code磕一个，但不幸的是不能找回插值 $Q_c$ 的数据就能很难完全再次做出和之前一样的 $Q_c$ 空间分布图了，并且剖面的插值也不能重现，我不知但我什么时候还能做出大差不差的结果了。
27记 （2026年7月6日 ~ 2026年7月12日） 这周可能是本学期的倒数第二周吧，教务说是要在16日左右对学生宿舍进行整修，所以就赶着我们回家。
- 很庆幸，目前为止那些被我大意删除的绘图脚本和中间数据基本得到了恢复和重制，除了上周多次尝试重现的 $Q_c$ 空间分布数据😥。在这期间，VS Code的自动备份文件机制和Google AI Studio的历时会话保存功能。
- 截至到撰写，我已经尝试了多种距离分析和空间连接（5km-20km、1km-2km）的克里金插值（参数目前保持默认）结果重现组合，但结果均不理想。我已经完全忘记了几个月前我所设置的参数了，并且当时没有对此有所文字记录，所以我只能一点点的对往日的结果摸索着重现。至于能否做出相同的结果，我不敢保证，但我希望我能尽快。
- 《地震地质》的两位外审专家，其中一位已经完成了他的审稿，一旦另一位也审稿结束，由于我还不能完成复刻 $Q_c$ 的空间分布图及其剖面图，要是让我大修包括数据在内的工作，不知道我该如何面对接下来的处理。

* 一些疑问

✨🗝据经验及地震的实际观测，为消除衰减的S波对尾波记录的影响，通常将尾波窗口起始时间取为2倍的S波走时的时间。

这里建议 将尾波窗口起始时间取为2倍的S波走时的时间 ，为什么姜和梁是以发震时间和P波到时为基础，计算理论2倍S波到时 $2\sqrt{3}(T_p-T_o)$ ，而不是手动拾取S波（或者使用深度学习的方法？）再取2倍S波走时的时间？这一近似的S波到时是否会对尾波衰减 $Q_c$ 的计算产生某些影响？

答复： 原文并没有在这一处理上进行详细解释说明。是因为人工拾取P波较拾取S波准确且容易得多吗？是的
✨🗝在姜的数据文件中，在台站信息表zdyinfonew、station中的B63D、B741都代表zdy07？

答复： （2025年10月20日） 经验证（邹锐的论文基于气枪主动源的祁连山地区地下介质衰减变化中的台站位置图示）：B63D是zdy07，B741是zdy04。

如果后续需要使用姜已处理的数据，对这一部分的数据需要额外注意。

（2025年12月29日） 依据往日的台站信息表格（主动源台站信息.xls），B63D是zdy04，B741是zdy07，但经过我之前使用arcfetch和mseed2sac程序提取得到的SAC文件文件名来看，与往日的表格正好相反，其余编号都与之一致。

(2026年1月4日） 发现我拿到的主动源台站编号混乱，多个台站的数据格式转换出来得到的数据头段变量有一定的偏差，在后续处理中可能还需要进行额外的文件名识别。
✨🗝姜和梁对主动源数据的波形截取方式不一样，姜在提取特定地震事件波形时就已经裁剪了固定时间长度的波形，而梁则是在后续写入事件信息时使用cut命令来进行裁剪波形。
- 姜提前对地震发震时间减去10秒（见疑问6），使用arcfetch程序来提取出相对发震时间 -10~190s 的波形；
- 还未找到梁数据文件夹中关于主动源数据提取的脚本或方式描述。
答复： （2025年11月19日） 对于波形的裁剪，我采用姜的处理方式。此种方式得到SAC文件名会比实际发震时间早10s，故在写入发震参考时间 o 时可以根据文件名的时间加上10s即可，无需再遍历地震目录的时间。

（2025年11月19日下午） 又尝试通过arcfetch利用发震时间直接截取300s，根据上午的测试，不提前10s一般也能完整的得到整个波形记录。
✨🗝如何对提取的波形数据进行处理？

答复： 目前已根据地震目录，将主动源的数据截取了从发震时刻 $o$ 截取地震发生前10s，及之后200s的地震波形：
- 编辑SAC头字段写入发震时刻 $o$ 和标记P波到时 $T_1$ ？
- 去均值、去线性趋势和波形尖灭？rmean;rtrend;taper
已知，地震尾波衰减研究是可以不需要进行去仪器响应的，而波形筛选的要求是？

没有定量的条件吗，仅凭眼看波形质量？也许是。

对于不合适的波形，请参照SAC_Docs手册的5.17 质量控制一节，比如在任意时间点标记t9并将该类标记的数据移动到专用文件夹saclst t9 f *.SAC | awk '$2>0 {print "mv", $1, "BAD/"}' | sh。
✨🗝在姜的数据文件夹中，地震目录的时间和使用实际截取的波形数据时间不一致。

以\ZDY1_3\40\这一数据文件夹为例，30km内的地震目录文件30.EQT中显示2016年7月15时16分36秒60毫秒有一次地震事件，但在用于读取的发震时间文件w302（年积日）中却是2016年197日8时36分50秒000毫秒，这一时间处理（同一地震事件时分秒的不同）是如何进行的？

答复： （2025年11月18日） 文件夹中的bjtime2worldtime.sh脚本对北京时间和世界时进行了转换！将北京时间减去28800秒即可得到世界时！但在实际处理时，姜额外减去了10秒，是为了在截取地震波形数据时得到地震发生前10s的波形数据。我将该UTC时间转换功能整合到了convert_events_DoY.sh脚本中，启用--utc便可实现世界时的同步转换！
姜在选择尾波窗长LTW（流逝时间窗Lapse Time Window）时，为什么选择20、30、40s作为研究对象？

Kopnichev和Gao等人在实验中证明，尾波窗小于100 秒内观察到的尾波中主要成分是单次散射，而大于100秒时观察到的尾波中主要成分是多次散射。据实验分析，尾波窗长最小应为20秒，在实际应用中，大多选择20 - 40s。这里涉及一本外文书籍《Seismic wave propagation and scattering in the heterogeneous earth : Second edition》，需要了解一下姜在文中描述的据实验分析，尾波窗长最小应为20秒，在实际应用中，大多选择20 - 40s是怎么得到的？

另外，朱新运的博士论文《衰减、场地响应等地震波传播相关信息综合研究》中也提到以信噪比高于根号2截断尾波作为最大可用尾波段，对于记录质量较高的地震波以此标准截断的尾波窗可能很长尾波窗太长不可避免的引入了多次散射作用对此尾波窗长不宜太长 ，但文中并未过多提到尾波窗长的参数选择方法。尾波流逝时间是影响地震波尾波衰减参数最大的因素，反映到实际计算中表现为值随尾波流逝时间增大而增大。朱新运在其研究中，将流逝时间严格限定为60s。
✨🗝对于一些截取出来的200s波形数据，在同一条波形中有两个甚至以上个数的地震记录，比如在zdy35的2022年第7和第8日就有大量的多地震波记录，这一情况下应该如何处理？

答复： （2025年12月24日下午） 目前对于该类型的数据（主要是zdy35及其之后的台站数据）处理办法是，如果两个地震时间相距不是特别近则优先拾取第一个地震P波而第二个地震P波则无法拾取；如果时间相距比较远，则优先拾取第二个地震P波，且第一个地震P波也可以进行拾取。
✨🗝对于我重构的用于计算 $Q_c$ 的python脚本，目前有两个需要注意的地方。
- 在计算 $SNR$ 时，往日直接采用P波到时的两倍，而在当前则是2倍S波到时（2倍根号3）
- 在计算 $Q_c$ 时，往日代码采用i-1进行索引，而当前则直接使用全部索引（不减1）
该疑问描述请见10记部分的内容。

答复： （2026年1月6日下午） 在计算 $SNR$ 信噪比最后5秒的python代码处应当采用当前则是2倍S波到时（2倍根号3）这一原则，师兄师姐他们在此处有一个错误，需要注意！正确地采用 2倍根号3的P波走时(即2倍S波走时) 来计算信噪比，会比错误地采用 _{2倍P波走时} 得到的符合筛选条件的记录少。这是因为错误地采用 _{2倍P波走时} 时，用于计算信噪比的波记录大概率并不在尾波时间窗之内，导致这样得到的 $SNR$ 偏大，即同等筛选条件下（ $SNR\ge 3$ ）的记录会偏多。

关于此处的计算索引，目前测试结果来看，并不会有任何错误或影响，保持当前索引方式。
✨🗝为什么要使用 $Q_c{−1}$ 而不是 $Q_c$ 进行插值？

物理机制的依据：能量衰减的线性可加性。在地震学中，品质因子 $Q$ 是一个无量纲参数，它本身并不直接代表能量损失的绝对量。真正代表介质“内摩擦”或“衰减系数”（即地震波每传播一个波长所损失的能量比例）的是它的倒数 $Q^{−1}$ 。从衰减理论公式 $Q_{total}^{−1}=Q_{intrinsic}^{−1}+Q_{scattering}^{−1}$ 可以看出，物理上的衰减效应是严格按照 $Q^{−1}$ 来进行线性叠加的。空间插值（无论哪种算法）的数学前提是变量在空间上具有线性连续性。对 $Q$ 插值在物理上是错误的，必须对 $Q^{−1}$ 插值才符合波的能量耗散积分原理。

数学分布的依据：压制非正态分布畸变。 $Q_c$ 值在自然界中往往跨度极大（从几十到几千），其空间分布通常呈严重的偏态（Log-normal distribution）。如果直接对 $Q_c$ 插值，高达2000的极值点会像“黑洞”一样拉高周围的所有预测值，掩盖真实的衰减异常。而转换成 $Q^{−1}$ 后，数值分布更加平稳、趋于正态分布，完美契合地质统计学的要求。

Sato, H., Fehler, M. C., & Maeda, T. (2012). Seismic wave propagation and scattering in the heterogeneous earth. Springer. （国际公认的地震波散射圣经，明确指出衰减属性的物理表达为 $Q^{−1}$ ）

赵连锋, 谢小碧, 何熹, 等. (2022). 地震 Lg 波衰减成像方法, 算法, 数据处理流程及应用. 地球与行星物理论评, 53(6), 24.（国内顶级的衰减层析成像综述，强调衰减系数 $Q^{−1}$ 是空间积分与反演的基础）
✨🗝为什么要对 $Q_c^{−1}$ 进行渔网均值化，而不是直接平均 $Q_c$ ？

避免“高值掩盖低值”的严重数学偏倚（Nugget Effect Bias）。地质学上，我们最关心的是断裂带、流体、破碎带，这些区域表现为高衰减（极低的 $Q_c$ 值，例如 $Q_c=100$ ）。而稳定的岩石表现为低衰减（极高的 $Q_c$ 值，例如 $Q_c=1000$ ）。
- 如果一个网格内有两条射线，一条测出Q_c=100，一条测出Q_c=1000。
  - 错误做法（算术平均 $Q_c$ ）： $(100+1000)/2=550$ 。结果：这个网格看起来是个中等偏稳定的区域。 断层的低 $Q_c$ 异常被彻底抹杀了！
  - 正确做法（平均 $Q_c^{−1}$ ，即物理上的调和平均）： $(1/100+1/1000)/2=0.0055$ 。 将其倒数还原，真实的平均 $Q_c\approx181$ 。结果：该网格成功保留了高衰减（低 $Q_c$ ）的物理特征！
- 因此，对 $Q_c$ 直接求平均，会严重抹杀掉构造活跃区（低 $Q$ 异常）的特征，只有对 $Q_c^{−1}$ 求平均，才能真实保留介质的破碎信息。
Bao, X., Sandvol, E., Ni, J., et al. (2011). High resolution regional seismic attenuation tomography in eastern Tibetan Plateau and adjacent regions. Geophysical Research Letters, 38(16).（青藏高原东缘衰减成像的经典文献，其反演和网格平滑都是基于 $Q^{-1}$ 或衰减系数展开的）。

Pei, S., Liu, J., Ma, H., et al. (2010). Dynamic Variation of S-Wave Q Value Beneath Sichuan Yunnan, China. Chinese Journal of Geophysics, 53(7), 1639-1652.（在处理多射线交叉网格时，采用衰减率 $Q^{-1}$ 进行平滑分配）。

** 如何拟合频率依赖关系

在梁和姜的论文中，这对应于“计算得到的尾波 $Q_c$ 随频率变化关系”部分（通常形如 $Q_c = Q_0 f^\eta$ ）。实现这一步的核心思想是将非线性的幂函数关系转化为线性的对数关系，然后利用最小二乘法进行回归。

i. 对数线性化

假设尾波 $Q_c$ 与频率 $f$ 满足幂函数关系：$ Q_c(f) = Q_0 \cdot f^\eta $

$Q_0$ ：表示频率为 1Hz 时的 $Q_c$ 值，反映了该区域整体的衰减水平（构造越活跃， $Q_0$ 通常越低）
$\eta$ (eta)：频率依赖指数，反映了地下介质非均匀性的程度（构造越活跃， $\eta$ 通常越高）

为了求解 $Q_0$ 和 $\eta$ ，在公式两边取自然对数（或常用对数）： $\ln(Q_c) = \ln(Q_0 \cdot f^\eta)$ $\ln(Q_c) = \ln(Q_0) + \eta \cdot \ln(f)$

这实际上就是一个一元一次线性方程 $Y = ax + b$ ：

$Y = \ln(Q_c)$
$X = \ln(f)$
斜率 $a = \eta$
截距 $b = \ln(Q_0)$

实施目标： 通过线性回归求出斜率和截距，进而得到 $\eta$ 和 $Q_0$ （ $Q_0 = e^b$ ）。

ii. 详细步骤

前提：已计算得到某台站（或某区域）在所有地震事件下、各个频率（1.5, 3, …, 24 Hz）的 $Q_c$ 值。

步骤 1：数据清洗与平均 (Data Aggregation)

由于单个地震计算出的 $Q_c$ 值可能存在波动，不能直接对成百上千个散点进行拟合（ $R^2$ 会很低且权重不好控制），而应该先计算每个频率点下的平均值和标准差。

将所有数据按中心频率分组（1.5Hz 组, 3Hz 组…）；
剔除每个组内的极端异常值（例如超过 3 倍标准差的数据）；
计算每个频率点下 $Q_c$ 的均值 (Mean) 和 标准差 (Std)。

步骤 2：对数变换 (Log Transformation)

准备两个数组：

$X_{log} = [\ln(1.5), \ln(3), \ln(6), \dots, \ln(24)]$
$Y_{log} = [\ln(Q_{c\_mean}^{1.5}), \ln(Q_{c\_mean}^{3}), \dots, \ln(Q_{c\_mean}^{24})]$

步骤 3：加权线性回归 (Weighted Linear Regression)

虽然普通最小二乘法（OLS）也可以，但考虑到不同频率下的数据质量不同（例如低频部分信噪比可能低，样本离散度大），进阶的做法是使用加权最小二乘法，权重为标准差的倒数。

拟合 $X_{log}$ 和 $Y_{log}$ ，得到：

斜率 $a$ $\rightarrow$ 即为 $\eta$
截距 $b$ $\rightarrow$ 计算 $Q_0 = e^b$

步骤 4：误差估计与绘图

计算拟合优度 $R^2$
绘制双对数坐标图：横轴为频率 $f$ ，纵轴为 $Q_c$
画出平均值点、误差棒（标准差）以及拟合的直线

iii. 代码实现

import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import r2_score

# ==========================================
# 1. 数据录入（每个频率下的观测值，二维列表/数组）
# ==========================================
# 每个子列表为一个频率下的所有观测值
freqs = np.array([1.5, 3, 6, 9, 12, 18, 24])
# 从文本文件读取观测值
def read_txt_data(filename):
    with open(filename, 'r', encoding='utf-8') as f:
        return [float(line.strip()) for line in f if line.strip()]

# 读取各个频率的观测值
time_window = 40  # seconds
qc_data = [
    read_txt_data(f'./{time_window}s/{time_window}s1p5Hz.txt'),  # 1.5Hz
    read_txt_data(f'./{time_window}s/{time_window}s3Hz.txt'),    # 3Hz
    read_txt_data(f'./{time_window}s/{time_window}s6Hz.txt'),    # 6Hz
    read_txt_data(f'./{time_window}s/{time_window}s9Hz.txt'),    # 9Hz
    read_txt_data(f'./{time_window}s/{time_window}s12Hz.txt'),   # 12Hz
    read_txt_data(f'./{time_window}s/{time_window}s18Hz.txt'),   # 18Hz
    read_txt_data(f'./{time_window}s/{time_window}s24Hz.txt')    # 24Hz
]
# 自动计算均值和标准差
# Q_c^{-1}：均值采用调和平均数 (Harmonic Mean)：先求倒数的平均，再倒数回来
# 标准差使用误差传递公式（Delta Method）: sigma_Q_c = sigma_inv * (Q_c)^2

# 先计算1/Q_c的统计量
qc_inv_data = [1.0 / np.array(arr) for arr in qc_data] # 获取某频率下所有的 Qc^-1
mu_inv = np.array([np.mean(arr) for arr in qc_inv_data]) # Step 1: 算术平均值 mu_inv
sigma_inv = np.array([np.std(arr, ddof=1) for arr in qc_inv_data]) # Step 1: 标准差 sigma_inv

# 通过调和平均计算Q_c的均值
qc_means = 1.0 / mu_inv # Step 2: 真实的区域平均 Qc = 1 / mu_inv

# 使用误差传递公式计算Q_c的标准差
qc_stds = sigma_inv * (qc_means ** 2) # Step 3: 根据公式 sigma_Qc = sigma_inv * (Qc)^2

print("观测数据点数:", [len(arr) for arr in qc_data])
print("观测数据均值:", qc_means)
print("观测数据标准差:", qc_stds)

# ==========================================
# 2. 数据预处理与权重计算
# ==========================================
# 转换到对数域
x_log = np.log(freqs)
y_log = np.log(qc_means)

# 计算对数域的标准差 (误差传播: d(lnx) = dx/x)
# sigma_log_y = sigma_y / y
log_stds = qc_stds / qc_means

# 计算权重 (Weights = 1 / sigma)
# 这一步保证了相对误差小的数据点（更可靠的点）在拟合中占比更重
weights = 1.0 / log_stds

# ==========================================
# 3. 加权线性拟合
# ==========================================
# 使用 numpy.polyfit 进行拟合
# deg=1 表示拟合直线 y = kx + b
# cov=True 返回协方差矩阵，用于计算拟合参数的误差
coeffs, cov_matrix = np.polyfit(x_log, y_log, deg=1, w=weights, cov=True)

eta_fit = coeffs[0]          # 斜率即为 eta
ln_q0 = coeffs[1]            # 截距为 ln(Q0)
Q0_fit = np.exp(ln_q0)       # 反解 Q0

# 计算误差 (Standard Errors)
perr = np.sqrt(np.diag(cov_matrix)) # 对角线开根号得到参数标准误
eta_error = perr[0]
ln_q0_error = perr[1]
# Q0 的误差传播: delta_Q0 = Q0 * delta_ln_Q0
Q0_error = Q0_fit * ln_q0_error 

# 计算拟合优度 R^2
y_pred_log = eta_fit * x_log + ln_q0
r2 = r2_score(y_log, y_pred_log, sample_weight=weights)

# ==========================================
# 4. 输出结果
# ==========================================
print("="*40)
print("拟合结果 Report")
print("="*40)
print(f"Q0  = {Q0_fit:.2f} ± {Q0_error:.2f}")
print(f"eta = {eta_fit:.4f} ± {eta_error:.4f}")
print(f"R^2 = {r2:.4f}")
print("-" * 40)
print(f"拟合公式: Qc(f) = ({Q0_fit:.1f}±{Q0_error:.1f}) * f^({eta_fit:.2f}±{eta_error:.2f})")
print("="*40)

# ==========================================
# 5. 专业绘图 (对数坐标系)
# ==========================================
# 设置字体：英文使用Times New Roman，中文使用宋体
plt.rcParams['font.serif'] = ['Times New Roman', 'SimSun']
plt.rcParams['font.family'] = 'serif'
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题
# 设置数学公式字体为Times New Roman风格
plt.rcParams['mathtext.fontset'] = 'stix'  # STIX字体与Times New Roman相似
plt.rcParams['mathtext.rm'] = 'Times New Roman'
plt.rcParams['mathtext.it'] = 'Times New Roman:italic'
plt.rcParams['mathtext.bf'] = 'Times New Roman:bold'

plt.figure(figsize=(7, 5), dpi=210)

# 绘制所有原始观测点（散点）
for i, f in enumerate(freqs):
    plt.scatter([f]*len(qc_data[i]), qc_data[i], color='#5DADE2', alpha=0.5, s=30)

# 绘制拟合曲线
f_smooth = np.linspace(1, 30, 100)
qc_fit_curve = Q0_fit * (f_smooth ** eta_fit)
plt.plot(f_smooth, qc_fit_curve, linestyle='--', linewidth=4, color='#FF5252',
         dash_capstyle='round',
         label=f'Fit: <span class="math math-inline" data-math-placeholder="33"></span> (<span class="math math-inline" data-math-placeholder="34"></span>)')

# 绘制均值和误差棒
plt.errorbar(freqs, qc_means, yerr=qc_stds, fmt='o', color='#34495E', 
             ecolor="#70DE73", elinewidth=1.5, capsize=4, 
             label='Mean ± Std')

# 设置坐标轴为对数刻度
plt.xscale('log')
plt.yscale('log')

# 设置刻度显示格式
from matplotlib.ticker import ScalarFormatter
plt.gca().xaxis.set_major_formatter(ScalarFormatter())
plt.xticks([1.5, 3, 6, 9, 12, 18, 24], ['1.5', '3', '6', '9', '12', '18', '24'])
# 手动设置Y轴刻度位置，显示更多标注
plt.yticks([10, 100, 1000, 10000])

# 标签与美化
plt.xlabel('Frequency (Hz)', fontsize=11)
plt.ylabel('<span class="math math-inline" data-math-placeholder="35"></span>', fontsize=11)
plt.title(f'Frequency Dependence of Coda <span class="math math-inline" data-math-placeholder="36"></span>\n with {time_window}s Time Window', fontsize=13)
plt.grid(True, which="both", ls="-.", alpha=0.3)
plt.legend(fontsize=10)

plt.show()

*** 如何使用克里金法进行空间插值

核心逻辑是：离散点 -> 空间聚合(降噪/平滑) -> 稀疏点 -> 全局插值。

i. 准备离散点

已通过使用数据管理工具 > 要素 > XY转线生成了地震射线路径

1. 打开工具箱：

Data Management Tools (数据管理工具) > Sampling (采样) > Generate Points Along Lines (沿线生成点)。

2. 设置参数：

Input Features (输入要素)：射线（线图层）。
Output Feature Class (输出要素类)：设置输出的点文件路径。
Point Placement (点位置)：选择 DISTANCE (按距离)。
Distance (距离)：这是最关键的参数。输入采样间隔，例如 5 Kilometers (5千米) 或 0.05 Decimal Degrees (取决于你的坐标系)。间隔越小，点越密，插值越精细，但计算量越大。
Include End Points (包含端点)：勾选。这能保证台站位置和震源位置都被采样到。

3. 检查坐标系：

确保 Ray_Points 使用的是投影坐标系（如 UTM），而不是地理坐标系（经纬度）。
原因：后续建立渔网需要以“米/千米”为单位，经纬度的度数在不同纬度代表的距离不同，会导致网格变形。
工具：Project (投影)。

ii. 消除射线路径

这是消除“条带效应”的关键步骤。

1. 创建渔网 (Create Fishnet)

工具箱：Data Management Tools > Sampling > Create Fishnet
关键参数设置：
- Output Feature Class：命名为 Grid_Mesh.shp。
- Template Extent (模板范围)：Ray_Points 图层（或研究区边界）。
- Cell Size Width/Height (像元大小)：核心参数！
  - 建议设置为 20km x 20km 或 0.2 度。
  - 经验法则：网格大小应略大于射线间的平均间隙。如果太小，还是会有空洞；如果太大，分辨率会降低。
- Geometry Type：选择 POLYGON。
结果：得到覆盖研究区的整齐方格网。

2. 空间连接 (Spatial Join)

工具箱：Analysis Tools > Overlay > Spatial Join
关键参数设置：
- Target Features (目标要素)：Grid_Mesh (刚才生成的网格)。
- Join Features (连接要素)：Ray_Points (密集的射线点)。
- Output Feature Class：命名为 Grid_With_Qc.shp。
- Join Operation：JOIN_ONE_TO_ONE。
- Field Map of Join Features (字段映射) —— 最重要的一步！：
  1. 在列表中找到 Qc 字段。
  2. 右键点击 Qc，选择 Merge Rule (合并规则)。
  3. 选择 Mean (平均值)。
  - (这一步的作用是：如果一个网格里有100个点，计算这100个点的Qc平均值赋给这个网格)。
- Match Option：CONTAINS (包含)。
结果：得到一个新的网格图层，属性表中有一个 Qc 字段（实际是平均值）。很多没射线的网格该字段为空（Null）。

3. 筛选有效网格 (Select & Export)

由于研究区外有很多空网格，需要剔除。
打开 Grid_With_Qc 属性表 -> 按属性选择 -> Qc IS NOT NULL。
右键图层 -> Data -> Export Features -> 保存为 Valid_Blocks.shp。

4. 面转点 (Feature To Point)

插值工具通常需要点作为输入。
工具箱：Data Management Tools > Features > Feature To Point
Input Features：Valid_Blocks.shp。
Output Feature Class：Block_Centers.shp。
结果：一组分布均匀、不仅保留了空间变化趋势，还消除了单条射线噪声的高质量采样点。

iii. 普通克里金插值法

现在使用处理过的 Block_Centers 进行插值。

1. 普通克里金插值 (Kriging)

工具箱：Spatial Analyst Tools > Interpolation > Kriging (或者用地统计向导)。
Input Point Features：Block_Centers。
Z value field：Qc。
Semivariogram Model (半变异函数模型)：
- 选择 Spherical (球状) 或 Gaussian (高斯)。
- 此时因为数据已经是平滑过的中心点，Nugget (块金值) 可以设得很小甚至为0。
Search Radius (搜索半径)：
- 设置为 Variable (可变)，点数设为 12 左右即可（因为现在的点已经是代表区域的平均值了，不需要搜太远）。
Output Cell Size：可以设置得比之前的网格小（例如 5km），以获得平滑的视觉效果。
Environments (环境设置)：
- Processing Extent：设置为研究区边界。
- Mask：设置研究区掩膜。