简介本资源面向地震学与深度学习交叉方向的学生及开发者提供一套基于Python实现的微震拾取模型完整工程可用于毕业设计、课程设计或项目开发。模型以16秒长、100Hz采样的三分量地震波形为输入经标准化后转为(1,3,1600)张量输出(1,2,1600)的P/S波初至概率无需滤波即可在2级以下地震事件中保持较好鲁棒性。压缩包共21个文件约1.57MB包含6个py源码文件模型定义、数据加载、训练与评估脚本、3个pt权重文件、2个ipynb实验笔记、8张png结果图及1份md说明文档覆盖从训练到可视化评估的完整流程。已有56人学习下载。读者可获取可直接运行的源码、预训练权重与项目文档理解ResUnet结构在微震拾取中的落地方式并在此基础上迁移至自身数据或延伸改进。1. 微震拾取模型到底在做什么从一段波形到一次触发井下微震监测的原始数据本质上是几十路检波器连续采样出来的时间序列。值班人员真正关心的不是波形好不好看而是「几点几分几秒在哪个位置发生了一次微震事件」。微震拾取模型要解决的就是这条链路里最靠前、也最耗人力的一环从连续波形里判断出 P 波初至时刻并把它和噪声、钻机干扰、电气脉冲区分开。传统做法靠 STA/LTA 这类短长时窗能量比算法参数调得好在信噪比高的台站上够用但换一个工作面、换一批传感器阈值就得重调误拾和漏拾同时存在。基于深度学习加 Python 实现的微震拾取模型思路是把这件事当成逐采样点的二分类或语义分割问题输入一段多通道波形输出每个采样点属于 P 波到时的概率再取概率峰值作为拾取结果。这套方案适合三类人做毕业设计需要完整可跑通链路的同学、课程设计想找一个有真实数据形态的题目、以及现场做监测系统开发、想把人工拾取替换成自动拾取的工程师。它不要求你有 GPU 集群一张消费级显卡甚至纯 CPU 都能把最小闭环跑起来关键在于数据组织和标签质量而不是模型有多深。2. 数据准备与标签构造微震拾取模型能不能用八成看这一步2.1 微震数据的三种常见来源与格式差异做这个模型第一件翻车的事往往不是模型而是数据。微震数据来源大致三类格式和处理方式差别很大。第一类是标准地震学格式比如 miniSEED、SAC。这类数据自带采样率、起始时间、台站信息用 ObsPy 读取最省事缺点是文件切分碎一个台站一天可能几百个文件需要先做合并和重采样。第二类是采集系统导出的自定义二进制或文本常见的是每通道一个文件或者一个文件里按通道交错存放。这类数据没有元信息采样率、通道顺序、字节序都得从采集配置里确认字节序搞反了波形会变成一片噪声这是新手最容易踩的坑。第三类是已经切好的事件段通常是 npy 或 mat 格式每个文件是一段固定长度的多通道波形附带一个到时标签。这类数据拿来训练最直接但要注意它可能已经被滤波或去均值过如果训练和推理的预处理不一致模型上线就会掉点。我一般会先写一个统一的数据探查脚本把采样率、通道数、时长、幅值范围打印出来确认没有异常通道再往下走。import obspy import numpy as np # 读取单个 miniSEED 文件确认元信息 st obspy.read(data/station01.mseed) print(通道数:, len(st)) print(采样率:, st[0].stats.sampling_rate) print(起始时间:, st[0].stats.starttime) print(时长(秒):, st[0].stats.endtime - st[0].stats.starttime) # 转成 numpy 数组形状为 (通道, 采样点) data np.stack([tr.data for tr in st]).astype(np.float32) print(数据形状:, data.shape, 幅值范围:, data.min(), data.max())这段代码的作用是把元信息和数值一次性看清楚。sampling_rate决定后面窗口长度怎么换算成秒data.shape决定模型输入是单通道还是多通道。如果幅值范围出现 1e10 这种量级说明字节序或数据类型解析错了先解决这个再谈训练。2.2 标签从哪来人工拾取、模板匹配与半自动标注微震拾取是逐点标注任务标签就是每个采样点是不是 P 波初至。真实项目里标签来源有三种。人工拾取最准但成本高一个熟练人员一天也就标几百条。模板匹配适合台站固定、事件波形相似的场景用一个已知事件的波形去和连续数据做互相关超过阈值的位置自动打标再人工复核。半自动标注是现在比较务实的做法先用 STA/LTA 或一个粗糙模型跑一遍把候选段切出来人工只做「确认或微调到时」效率能提升好几倍。标签格式建议统一成「事件段 到时索引」。比如一段 6000 采样点的波形P 波到时在第 1800 点标签就是 1800。如果做逐点分类就把 1800 附近若干点标为 1其余为 0如果做回归就直接预测这个索引。分类对类别不平衡更鲁棒回归输出更直接我一般先用分类跑通再考虑回归。2.3 把连续波形切成训练样本窗口长度与正负样本比例模型输入是固定长度窗口窗口长度要覆盖一个完整 P 波加一段背景噪声。经验值是 2 到 4 秒采样率 500 Hz 的话就是 1000 到 2000 点。太短会切掉 P 波尾部太长会引入无关噪声、增加计算量。正负样本比例是另一个关键。真实数据里微震事件稀疏如果随机切窗口负样本可能占 95% 以上。直接训练模型会倾向于全预测为负。常见做法是正样本窗口以到时为中心切负样本窗口从无事件段随机切把正负比控制在 1:1 到 1:3 之间。也可以在损失函数里用类别权重但样本层面先平衡更直观。import numpy as np def make_windows(data, pick_idx, win_len2000, pos_ratio2): data: (通道, 采样点); pick_idx: P波到时索引列表 n_ch, n_pt data.shape half win_len // 2 pos, neg [], [] for idx in pick_idx: if idx - half 0 or idx half n_pt: continue pos.append(data[:, idx - half: idx half]) # 负样本从远离所有到时的位置随机切 pick_set set(pick_idx) tries 0 while len(neg) len(pos) * pos_ratio and tries 100000: tries 1 start np.random.randint(0, n_pt - win_len) center start half if all(abs(center - p) win_len for p in pick_set): neg.append(data[:, start: start win_len]) X np.stack(pos neg).astype(np.float32) y np.array([1] * len(pos) [0] * len(neg), dtypenp.float32) return X, ywin_len是窗口采样点数要和采样率一起换算成秒来确认合理性。pos_ratio控制负样本相对正样本的倍数2 表示负样本是正样本两倍。abs(center - p) win_len保证负样本窗口和任何事件都不重叠避免把事件尾部当噪声。这段逻辑跑完X的形状是 (样本数, 通道数, 窗口长度)可以直接喂给 PyTorch 或 TensorFlow。3. 模型选型与训练从 1D-CNN 到 U-Net 的取舍3.1 为什么微震拾取常用 1D-CNN 而不是直接上 Transformer微震波形是典型的一维时序信号局部波形特征P 波起跳的陡变、频率成分变化比长程依赖更重要。1D-CNN 的感受野通过堆叠卷积层就能覆盖几百个采样点参数量小、训练快、对数据量要求低非常适合毕业设计和中小规模项目。Transformer 在长序列建模上有优势但微震拾取里序列长度动辄几千点自注意力计算量是平方级而且需要大量数据才能训好。除非你有几十万条标注事件否则 1D-CNN 或 U-Net 结构的性价比更高。常见做法是编码器用几层卷积加池化解码器用上采样恢复分辨率最后逐点输出概率这就是 PhaseNet 那一类结构的核心思路。选型时还要考虑输入通道数。单台三分量用 3 通道输入台阵多通道用 N 通道输入第一层卷积的in_channels要对应改。通道数不是越多越好通道间不同步反而会干扰先确认各通道时间对齐再做多通道融合。3.2 一个能跑通的 1D-CNN 拾取网络与训练循环下面是一个最小可用的逐点分类网络输入 (batch, 通道, 窗口长度)输出每个采样点的概率。import torch import torch.nn as nn class PickNet(nn.Module): def __init__(self, in_ch3, base16): super().__init__() self.encoder nn.Sequential( nn.Conv1d(in_ch, base, 7, padding3), nn.ReLU(), nn.MaxPool1d(2), nn.Conv1d(base, base * 2, 5, padding2), nn.ReLU(), nn.MaxPool1d(2), nn.Conv1d(base * 2, base * 4, 3, padding1), nn.ReLU(), ) self.decoder nn.Sequential( nn.Upsample(scale_factor2, modenearest), nn.Conv1d(base * 4, base * 2, 3, padding1), nn.ReLU(), nn.Upsample(scale_factor2, modenearest), nn.Conv1d(base * 2, base, 3, padding1), nn.ReLU(), nn.Conv1d(base, 1, 1), ) def forward(self, x): return self.decoder(self.encoder(x)).squeeze(1)in_ch要和数据通道数一致base是基础通道数显存不够就调小。编码器两次池化把长度缩到四分之一解码器再上采样回来保证输出和输入等长。最后一层Conv1d(base, 1, 1)把特征压成单通道 logits配合BCEWithLogitsLoss使用。训练循环里要注意两点一是输入要做归一化按窗口减均值除标准差避免幅值差异导致梯度不稳二是正样本点远少于负样本点损失函数要加pos_weight。from torch.utils.data import DataLoader, TensorDataset def normalize(X): mean X.mean(axis2, keepdimsTrue) std X.std(axis2, keepdimsTrue) 1e-6 return (X - mean) / std # X: (N, C, L) 波形, y: (N, L) 逐点标签 X normalize(X) ds TensorDataset(torch.tensor(X), torch.tensor(y)) dl DataLoader(ds, batch_size32, shuffleTrue) model PickNet(in_chX.shape[1]) opt torch.optim.Adam(model.parameters(), lr1e-3) pos_weight torch.tensor([20.0]) # 正样本少加权 loss_fn nn.BCEWithLogitsLoss(pos_weightpos_weight) for epoch in range(30): model.train() total 0 for xb, yb in dl: opt.zero_grad() logits model(xb) loss loss_fn(logits, yb) loss.backward() opt.step() total loss.item() print(fepoch {epoch}, loss {total / len(dl):.4f})pos_weight是最需要调的参数。正样本占比 5% 左右时20 是个合理起点如果模型输出全是 0就继续加大如果到处误拾就减小。lr用 1e-3 起步loss 震荡就降到 1e-4。训练轮数不用太多这类任务通常 20 到 50 轮就收敛过拟合了看验证集 loss 早停。3.3 训练完怎么判断模型真的学到了 P 波光看 loss 下降不够要拿验证集做逐事件评估。对每个事件段取模型输出概率最大的位置作为拾取点和人工标签比误差。误差在 10 个采样点以内500 Hz 下约 20 毫秒就算合格。如果误差分布是双峰说明模型在某些事件上系统性偏移通常是标签本身不一致或预处理不统一。还要看误拾率在纯噪声段上跑模型统计有多少窗口输出了高概率峰值。误拾率高说明负样本不够或pos_weight太大。这两个指标一起看才能判断模型能不能上现场。4. 推理部署与现场落地从脚本到能用的拾取服务4.1 连续数据流式推理滑窗、重叠与去重训练是切好的窗口现场是连续数据流必须做滑窗推理。窗口按固定步长向前滑动步长小于窗口长度以保证重叠避免事件正好落在窗口边界被切掉。每个窗口输出一条概率曲线重叠部分取平均或取最大最后在整条概率曲线上做峰值检测。def stream_pick(model, data, win_len2000, step500, thr0.5): data: (通道, 总采样点)返回拾取点索引列表 model.eval() n_ch, n_pt data.shape prob np.zeros(n_pt) count np.zeros(n_pt) with torch.no_grad(): for start in range(0, n_pt - win_len 1, step): seg data[:, start: start win_len] seg (seg - seg.mean(axis1, keepdimsTrue)) / (seg.std(axis1, keepdimsTrue) 1e-6) x torch.tensor(seg[None]).float() p torch.sigmoid(model(x)).numpy()[0] prob[start: start win_len] p count[start: start win_len] 1 prob prob / np.maximum(count, 1) picks [] for i in range(1, n_pt - 1): if prob[i] thr and prob[i] prob[i-1] and prob[i] prob[i1]: picks.append(i) return picks, probstep是滑动步长500 表示每次前进 500 点重叠 1500 点重叠越多越稳但越慢。thr是峰值阈值现场宁可先调低再人工复核也不要漏掉事件。prob数组建议存下来方便事后回看模型在哪些位置犹豫这是调参和排查的主要依据。4.2 阈值、后处理与误拾抑制的工程参数模型输出的是概率真正决定拾取结果的是后处理参数。除了峰值阈值还要加两个约束最小峰间距和最小持续时间。最小峰间距防止一个事件被拾成多个一般设 0.5 秒对应采样点数最小持续时间要求概率超过阈值连续若干点滤掉单点毛刺。参数典型值作用调大后果调小后果峰值阈值0.3~0.5判定是否为拾取漏拾增多误拾增多最小峰间距0.5 秒抑制重复拾取密集事件被合并单事件多次拾取最小持续点数5~10滤除毛刺弱事件被滤掉噪声被当事件滑窗步长窗口 1/4控制重叠计算变慢边界事件丢失这些参数没有万能值要拿现场数据跑一批统计误拾和漏拾再折中。我一般会保留一份「模型概率 人工标签」的对照表每次调参都回看这张表避免凭感觉改。4.3 把模型接进现有监测系统的两种方式落地方式取决于现有系统。一种是离线批处理定时把采集系统导出的数据拉到服务器跑一遍推理把拾取结果写进数据库值班界面读库展示。这种方式改动小适合先验证效果。另一种是在线服务模型封装成 HTTP 或 gRPC 接口采集端实时推数据服务端返回拾取结果。这种方式延迟低但对稳定性和资源占用要求高要处理断流、乱序、重复推送。常见做法是先用离线方式跑一两个月确认误拾率可接受再考虑上线。不管哪种方式都要保留原始波形和模型概率出问题时能回放。微震拾取这种事没有后悔药只有可追溯的日志。5. 避坑与排查微震拾取模型最常见的五类翻车5.1 训练 loss 一直不降模型输出全是背景现象训练几十轮loss 卡在某个值不动推理时概率曲线几乎全为 0。原因正样本占比过低pos_weight没设或设得太小模型学到「全预测负」就能拿到很低的 loss。解决先统计训练集正负点比例把pos_weight设成负正比的量级比如正样本占 3% 就设 30 左右。同时检查标签是否真的对齐标签整体偏移几十点也会导致模型学不到。5.2 验证集误差很小现场误拾却很多现象验证集上误差十几毫秒拿到现场数据一跑噪声段到处是拾取。原因训练数据的负样本和现场噪声分布不一致。训练时负样本多来自同一批台站现场可能换了传感器或增加了新的干扰源。解决从现场数据里切一批纯噪声段加入训练集做负样本重新微调。也可以提高峰值阈值和最小持续点数先压住误拾再逐步补数据。5.3 多通道输入时模型效果反而变差现象单通道训练效果尚可改成三分量或多通道后误差变大。原因通道间时间没对齐或者某个通道质量差、噪声大把有用信号淹没了。解决先做通道间互相关确认没有系统性时延再逐通道评估把明显异常的通道剔除或降权。多通道不是简单堆叠质量比数量重要。5.4 换采样率后模型完全失效现象500 Hz 数据训的模型直接用在 100 Hz 数据上拾取全乱。原因模型学到的波形特征和采样率强相关窗口长度对应的物理时间变了频率成分也变了。解决统一重采样到训练时的采样率再推理。如果必须用不同采样率就重新训练或做迁移微调不要指望一个模型通吃。5.5 推理速度跟不上实时数据现象离线跑没问题在线接入后延迟越积越多。原因滑窗重叠太大、模型太大、或者用了 CPU 推理。解决先减小重叠步长把步长从窗口的 1/8 调到 1/4再考虑模型剪枝或量化有 GPU 就用 GPU批处理多个窗口一起推理。实时场景下延迟和精度要一起权衡不能只盯精度。6. 进阶技巧用概率曲线反推模型到底在看哪里模型跑通之后真正拉开差距的是对概率曲线的解读。我习惯把模型输出的逐点概率和原始波形叠在一起看重点看三件事峰值位置是否稳定、峰值宽度是否合理、峰值前后有没有次峰。峰值位置稳定说明模型对这类事件有把握峰值宽度太窄可能是过拟合到某个尖锐特征换个台站就失效峰值前后有次峰说明模型在几个候选位置之间犹豫这时候人工复核能发现标签本身有歧义。更进一步可以做一个简单的扰动测试把输入窗口整体平移几十个采样点看峰值是否跟着平移。如果峰值不动说明模型学的是窗口内的绝对位置而不是波形特征这种模型换数据必翻车。这个测试不用改代码几行脚本就能跑。def shift_test(model, seg, shifts(-50, -20, 0, 20, 50)): seg: (通道, 窗口长度)观察峰值位置随平移的变化 model.eval() base_peak None for s in shifts: x np.roll(seg, s, axis1) x (x - x.mean(axis1, keepdimsTrue)) / (x.std(axis1, keepdimsTrue) 1e-6) with torch.no_grad(): p torch.sigmoid(model(torch.tensor(x[None]).float())).numpy()[0] peak int(np.argmax(p)) if base_peak is None: base_peak peak print(f平移 {s:4d} 点, 峰值位置 {peak}, 相对基准偏移 {peak - base_peak})如果相对偏移和输入平移量基本一致说明模型跟踪的是波形本身如果偏移量远小于平移量甚至不变就要警惕。这个测试花不了几分钟但能提前暴露很多现场问题。我自己做这类项目最大的习惯是任何一次调参或换模型都先把概率曲线存下来和上一版对比。模型指标好看不代表现场好用概率曲线才是那个不会骗人的黑匣子记录。希望帮到你。本文还有配套的精品资源点击获取