简介这套基于深度学习与Python的微震拾取模型资料包面向毕业设计、课程设计与项目开发场景适合需要快速搭建地震P波、S波初至自动拾取方案的读者。模型采用ResUnet结构输入16秒100Hz三分量波形并标准化为(1,3,1600)张量输出(1,2,1600)概率序列对噪声有一定鲁棒性无需滤波预处理也能在2级以下微震事件中取得良好识别效果。压缩包共21个文件大小1.57MB包含6个Python源码数据加载、模型定义、训练、评估等模块、3个预训练pt权重、8张可视化结果PNG、2个Jupyter演示/测试脚本及1份Markdown项目文档目录结构清晰便于按模块检索与二次开发。已有56人学习浏览整套代码经过严格测试可直接参照文档复现实验理解从数据预处理、模型训练到结果评估的完整链路并在此基础上按需调整网络结构或数据接口适合作为毕业设计、课程设计的实用参考。1. 微震波形初至拾取为什么值得换一种思路微震P波和S波初至时刻是震源定位、震级评估和破裂过程反演的第一道门槛。传统STA/LTA、AIC在信噪比偏低时误拾率快速上升对三分量数据还要手动核对效率很难提上去。这个基于深度学习python实现的微震拾取模型把问题重定义为“给一段16秒、100Hz的三分量波形输出每个采样点是P波或S波初至的概率”用ResUnet直接把长度为1600的时间序列映射成两条概率曲线。源码里包含训练脚本、评估脚本、模型定义和项目文档适合作为毕业设计、课程设计也适合生产中做预研。如果你在一线处理微震数据会关心特征怎么对齐、标签怎么构造、训练怎么收敛如果你是学生更关心整个pipeline能否立刻复现——这两件事下文都会落到具体代码。2. 从原始波形到训练样本100Hz三分量窗口的构建与标准化2.1 为什么固定16秒、100Hz和1600个点输入shape (1,3,1600) 不是随意定的。16秒窗口能覆盖绝大多数微震事件的P波和S波到时差100Hz是微震台网常见的采样率两者相乘得到1600个时间点。窗口太短S波可能落在窗口外窗口太长模型需要更多下采样层处理长程依赖训练成本上升。如果你的原始数据是其他采样率读取后需要先统一到100Hz再做截断。常见做法是先用 scipy.signal.resample 把波形重采样到100Hz然后围绕事件触发时刻截取16秒窗口from scipy.signal import resample def build_window(wave, orig_rate, target_len1600): # wave: (3, n_samples) 三分量原始波形 resampled resample(wave, int(wave.shape[1] * 100.0 / orig_rate), axis1) n resampled.shape[1] if n target_len: # 截取中间区域尽量保留事件附近波形 start (n - target_len) // 2 resampled resampled[:, start:start target_len] else: # 不足16秒则在末尾补零 pad target_len - n resampled np.pad(resampled, ((0, 0), (0, pad))) return resampled.astype(np.float32)这段代码先把任意采样率换算到100Hz再按目标长度截取或补零。注意axis1是按时间轴重采样三分量共享相同的重采样比例。start (n - target_len) // 2取中间窗口为的是在事件相对位置不确定时尽量保留初至前后的完整波形。实际操作中我不会在每个batch都执行重采样而是把原始数据预处理成独立的1600点npy文件之后dataloader只需要简单读取。2.2 三个分量独立标准化模型要求输入已经完成标准化。这里标准化的是“每个样本的每个分量”而不是全局归一化。原因很直接不同台站的仪器增益、背景噪声水平差异巨大全局统计会被少数强噪声台站拉偏按样本分量独立标准化后输入只剩下波形形态信息模型更容易学习P/S波的到时特征。def normalize_wave(wave): # wave: (3, 1600) 三分量波形 wave wave.astype(np.float32) mean wave.mean(axis1, keepdimsTrue) std wave.std(axis1, keepdimsTrue) 1e-8 return (wave - mean) / std1e-8是为了防止静默台站出现零方差导致除零。标准化参数必须在训练和推理时保持一致不能用训练集全局统计替代样本内统计否则模型输出概率会出现整体偏移。项目摘要里强调“对三分量做标准化作为输入”指的就是这种逐样本逐分量的z-score处理。2.3 标签用高斯峰而不是0/1尖峰模型输出是(1,2,1600)每个时间点代表属于P波或S波初至的概率但训练标签不能做成只在初至时刻等于1、其余全0的one-hot。原因是0.02秒的采样间隔本身有限真实人工标注也有几个采样点的抖动one-hot会让模型把“差几个点”和“差很多点”视为同等级错误收敛很困难。我通常把标签写成以初至时刻为中心的高斯峰def gaussian_label(p_pick, s_pick, length1600, sigma3.0): label np.zeros((2, length), dtypenp.float32) idx np.arange(length) for i, pick in enumerate([p_pick, s_pick]): if pick is None or pick 0: continue label[i] np.exp(-0.5 * ((idx - pick) / sigma) ** 2) return labelsigma3表示允许初至附近约0.03秒的灰度容差当标注精度在20ms以内时可以把sigma缩到2。这里的标签构造策略直接影响模型输出概率峰的宽度。如果直接使用one-hot常见的现象是训练损失下降很慢、输出概率几乎全为0因为模型很难从1600个位置里精确预测唯一一个正样本换成高斯峰之后初至附近的点都有正监督信号输出自然变成连续概率曲线。2.4 Dataloader与数据划分读取波形和标签的逻辑集中在dataloader.py。一个典型的PyTorch Dataset结构如下from torch.utils.data import Dataset class MicroSeismicDataset(Dataset): def __init__(self, index_file): self.df pd.read_csv(index_file) def __getitem__(self, idx): item self.df.iloc[idx] wave np.load(item[wave_path]) # (3, 1600) wave normalize_wave(wave) label gaussian_label(item[p_pick], item[s_pick]) return torch.from_numpy(wave), torch.from_numpy(label) def __len__(self): return len(self.df)配套的索引文件至少需要这几列列名类型说明wave_pathstr预处理后的npy路径shape(3,1600)p_pickfloatP波初至所在采样点可为nans_pickfloatS波初至所在采样点可为nansplitstr数据划分标记train/valid/test这里把split列放进索引文件因为实际项目切分通常按事件或台站进行避免同一个事件的不同台站同时出现在训练集和验证集中造成指标虚高。如果某个事件没有S波s_pick填nan对应通道标签全0即可。项目里的dataloader.py只要能正确处理nan整个数据准备阶段就基本结束。3. ResUnet的编解码结构与P/S双通道输出设计3.1 为什么用ResUnet而不是普通UNet微震波形是一维时间序列直接套2D UNet需要把波形reshape成二维矩阵破坏时间连续性。常见做法是把所有Conv2d改成Conv1d保留UNet的编解码结构。ResUnet在UNet基础上为每个卷积块加了残差连接这对初至拾取很有价值P波初至前后波形幅值变化往往只有几个采样点深层网络需要同时保留原始波形细节跳跃连接负责恢复分辨率残差连接则帮助梯度跨层流动。项目模型代码基于ResUnet实现也就是在编码器和解码器里都埋入残差块。3.2 编码器残差块与时间长度减半编码器第一个Stage保持1600个时间点之后每个Stage把时间长度减半通道数从3逐步扩张。残差块的标准写法如下class ResBlock1d(nn.Module): def __init__(self, in_ch, out_ch, stride1): super().__init__() self.conv1 nn.Conv1d(in_ch, out_ch, kernel_size3, stridestride, padding1) self.bn1 nn.BatchNorm1d(out_ch) self.conv2 nn.Conv1d(out_ch, out_ch, kernel_size3, stride1, padding1) self.bn2 nn.BatchNorm1d(out_ch) self.relu nn.ReLU(inplaceTrue) self.shortcut nn.Sequential() if stride ! 1 or in_ch ! out_ch: self.shortcut nn.Sequential( nn.Conv1d(in_ch, out_ch, kernel_size1, stridestride), nn.BatchNorm1d(out_ch) ) def forward(self, x): out self.relu(self.bn1(self.conv1(x))) out self.bn2(self.conv2(out)) return self.relu(out self.shortcut(x))stride2的卷积同时完成降采样所以不需要单独加池化层。残差连接里的shortcut用1x1卷积调整通道数和时间长度保证相加时shape一致。如果对初至高频细节敏感可以把第一个块的kernel_size从3换成5但后续Stage建议仍用3避免过强的边缘效应。我一般把中间Stage的卷积核保持3因为这个模型输入长度只有1600不需要很大的全局感受野。3.3 解码器与跳跃连接时间分辨率的恢复解码器负责把100个时间点的特征逐步上采样回1600点。上采样通常用转置卷积或线性插值项目采用转置卷积的常见实现既能学习上采样规律又能减少人工设计的插值核。结构上每一层解码器先上采样再把编码器对应Stage的输出拼过来Stage输入通道输出通道时间长度Encoder13321600Encoder23264800Encoder364128400Encoder4128256200Encoder5256512100Decoder1512256200Decoder2256128400Decoder312864800Decoder464321600跳跃连接把第2层的细节特征直接送到解码器倒数第二层让网络在恢复时间分辨率时不必完全依赖上采样产生的插值信息。注意表中Decoder通道数写的是拼接后的输入通道数实际代码里需要先拼接再卷积否则通道数不匹配。3.4 输出层两个通道的Sigmoid概率模型的最后是一个1x1卷积把32通道压缩成2通道后面接Sigmoid得到每个时间点属于P波和S波的概率。这里选择Sigmoid而不是Softmax因为P波和S波初至不是互斥事件同一个时间点既可能既不是P也不是S两个通道应该独立输出概率。class ResUnetMicroSeismic(nn.Module): def __init__(self, in_channels3, out_channels2, base_channels32): super().__init__() self.enc1 self._make_stage(in_channels, base_channels, stride1) self.enc2 self._make_stage(base_channels, base_channels * 2, stride2) self.enc3 self._make_stage(base_channels * 2, base_channels * 4, stride2) self.enc4 self._make_stage(base_channels * 4, base_channels * 8, stride2) self.enc5 self._make_stage(base_channels * 8, base_channels * 16, stride2) self.out_conv nn.Conv1d(base_channels, out_channels, kernel_size1) self.sigmoid nn.Sigmoid() def forward(self, x): # x: (B, 3, 1600) e1 self.enc1(x) # ... 编码器链 logits self.out_conv(decoder_out) # (B, 2, 1600) return self.sigmoid(logits)这段伪代码展示了整体骨架。base_channels32是通道数起点数据量很小可以降到16数据充足可以升到48或64。模型输出长度应当与输入完全一致如果训练时发现输出长度不是1600优先检查转置卷积的stride、padding和编码器下采样次数是否匹配。4. 训练策略与损失函数让模型学会数准初至点4.1 损失函数BCEWithLogits而不是CrossEntropy模型输出2个通道的独立概率训练时使用二分类交叉熵BCE更自然。P波和S波初至在时间上是独立的同一个时间点理论上不可能同时是P和S但两个通道之间并不存在互斥约束。如果用CrossEntropy模型会被迫让两个通道的概率加和等于1这和真实物理过程不符。项目里推荐使用BCEWithLogitsLoss内部把Sigmoid和BCE融合在一起数值上比分开算更稳定。import torch.nn as nn criterion nn.BCEWithLogitsLoss(pos_weighttorch.tensor([1.5, 2.5]))pos_weight用于处理正负样本不平衡一段16秒波形里初至附近的高斯标签显著为1的时间点通常只有几十个其余都是背景负类。两个通道的正样本比例不同S波初至更难拾取所以我给S通道更高的权重一般P通道取1.02.0S通道取2.04.0。如果模型的漏检率高于误检率可以适当增大pos_weight如果输出概率普遍偏高、出现大量假峰则调小。4.2 优化器与学习率策略用AdamW作为默认优化器它对权重衰减的处理比Adam更规范。初始学习率1e-3、批次大小32在大多数微震数据集上都能在50个epoch内看到收敛趋势。optimizer torch.optim.AdamW(model.parameters(), lr1e-3, weight_decay1e-4) scheduler torch.optim.lr_scheduler.ReduceLROnPlateau( optimizer, modemin, factor0.5, patience8 )ReduceLROnPlateau会在验证损失连续8个epoch不下降时把学习率减半。比固定步长衰减更省心适合初至拾取这种收敛曲线不平滑的任务。如果数据量很大也可以在训练前10个epoch做线性热身从很小的学习率逐步升到1e-3减少前期震荡。4.3 训练脚本里的关键参数项目里的configs.py集中放训练参数阅读和修改都方便。我一般会这样设置参数默认值说明epochs150数据少时可缩短到80batch_size32显存不足时降至16lr1e-3数据量大可提到2e-3weight_decay1e-4防止过拟合sigma3.0高斯标签宽度pos_weight_p1.5P通道正样本权重pos_weight_s2.5S通道正样本权重sigma在configs里出现意味着标签生成方式与训练参数同源。修改sigma后需要重新生成所有训练标签不能只改训练循环里的参数。验证指标除了损失还应该记录P波和S波分别的拾取误差中位数、拾取率。如果验证集S波拾取率偏低先提高pos_weight_s而不是盲目加深网络。4.4 训练循环与早停判断实际运行时命令通常只有一行python train.py --epochs 120 --batch-size 32 --device cudatrain.py内部会调用configs.py中的默认配置并保存每个epoch结束后在验证集上表现最好的模型到result/models/。早停判断不建议只看损失因为损失受到背景负样本主导可能出现损失下降但拾取率变差的情况。更可靠的做法是每个epoch结束后在验证集上跑一次快速评估统计绝对误差小于0.1秒的比例保存这个比例最高的模型。训练结束后的result/models目录里通常会有best_model.pt和last_model.pt评估和可视化脚本默认读取best_model。如果训练过程中发现验证损失出现明显回弹优先检查学习率是否过高或pos_weight设置是否极端。5. 概率曲线到初至时间评估脚本与工程落地技巧5.1 用visualize_and_evaluate.py复现评估结果项目提供了visualize_and_evaluate.py把模型输出直接画成波形和概率曲线的对照图。常见的运行方式是python visualize_and_evaluate.py \ --model result/models/best_model.pt \ --index-file data/test.csv \ --save-dir result/figs这个脚本会遍历测试集对每个样本输出三分量波形、P波概率曲线、S波概率曲线以及人工标注位置。看曲线的时候重点观察P波初至附近是否出现尖锐的单一峰S波概率曲线是否被P波响应污染。如果S波概率峰总出现在P波概率峰同一位置说明网络没有学到两者的到时差关系可以尝试在标签构造阶段增大S波通道权重或者增加S波通道高斯峰的sigma。5.2 从概率曲线提取初至时间模型直接输出的概率曲线不能当成初至时间还需要做峰值检测和平滑from scipy.ndimage import gaussian_filter1d from scipy.signal import find_peaks def extract_pick(prob, threshold0.35, min_distance5): prob gaussian_filter1d(prob, sigma1.5) peaks, props find_peaks(prob, heightthreshold, distancemin_distance) if len(peaks) 0: return None best peaks[np.argmax(props[peak_heights])] return bestthreshold0.35表示概率峰值至少达到0.35才认为是有效初至对低信噪比事件可以降到0.25但会增加误触发。distance5防止同一个初至被拆成相邻两个峰。平滑sigma选择1.5既能压住高频毛刺又不会把真正的初至峰拉平太多。如果处理的是强噪声事件我通常先带通滤波到240Hz再做峰值检测但模型训练时并没有使用滤波器所以滤波与否应该单独验证不要默认进行。5.3 滑动窗口推理处理长波形实际微震记录可能长达几分钟甚至更长不能直接塞进模型。这时把长波形切成16秒窗口相邻窗口重叠2秒逐窗推理后拼接概率曲线。窗口边缘的初至会被不完整截断概率通常很低建议丢弃每个窗口开头和结尾各0.5秒的预测。如果同一个初至出现在两个相邻窗口的重叠区保留置信度更高的峰值即可。当事件刚好落在窗口边界时模型可能给出一个很宽的平缓概率肩而不是明显的单峰遇到这种情况把窗口向后平移2秒重新预测往往能恢复尖锐的P波峰值。5.4 导出ONNX接入实时处理训练完成后把模型导出为ONNX方便用C或Java调用。导出时注意让模型处于eval模式输入固定为(1,3,1600)python -c import torch; mtorch.load(result/models/best_model.pt); m.eval(); dtorch.randn(1,3,1600); torch.onnx.export(m, d, pickers.onnx, input_names[wave], output_names[prob], opset_version13)这段命令把PyTorch模型导出为ONNX输入命名wave输出命名prob。导出后可以用onnxruntime加载半精度推理在CPU上也能跑得很快。需要注意的是导出时的输入尺寸必须是(1,3,1600)不要在导出时设置动态时间轴否则后续部署要额外处理序列长度变化。生成后的pickers.onnx可以直接接入台网数据流每16秒滑窗调用一次即可。本文还有配套的精品资源点击获取
