简介：这份资源面向具备一定MATLAB基础、希望将深度学习应用于地震信号分析的学生与研究人员，提供了一套用卷积神经网络完成地震等级预测的完整代码实现。包内共33个文件，以30个m脚本为核心，辅以2个mat数据文件和1个xlsx表格，压缩包约370KB，涵盖网络搭建、前向与反向传播、梯度校验、数据预处理及超参数组合搜索等模块，结构紧凑便于逐层阅读。已有351人学习下载，说明其在同类小众方向中具备一定参考价值。读者可借此理解CNN从滤波器初始化、卷积池化到全连接分类的完整链路，掌握地震时间序列的标准化、降噪与时频特征提取思路，并借助数值梯度检查与随机种子设定保证实验可复现，进而迁移到自身的地震风险评估或信号分类任务中。
1. 地震等级预测为什么值得用 CNN 做一遍
地震等级预测在工程地震学和台网日常运维里是个老问题。传统做法要么走 Gutenberg-Richter 经验关系，要么用震级-频度统计外推，能给出一个大致范围，但对单台站、短窗口、早期波形这类信息量有限的场景，误差经常大到没法直接用。近几年台网密度上来了，三分量波形采样率普遍到 100 Hz 甚至更高，单事件可用的时序点数从几百涨到几千，这就给数据驱动方法留出了空间。卷积神经网络 CNN 最擅长的恰好是从这种高维时序里自动抽局部模式——P 波初动、S 波到时、尾波衰减这些特征，人工设计滤波器要反复调，CNN 靠卷积核自己学。MATLAB 在这条链路上有天然优势：信号处理工具箱、深度学习工具箱、Wavelet 工具箱都在一个环境里，读 miniSEED、做带通滤波、切窗、归一化、搭网络、训练、导出，不用在 Python 和 C 之间来回倒。这篇面向的是手里有台网波形、想用 MATLAB 把 CNN 地震等级预测跑通的从业者，从数据准备一路讲到调参和踩坑，新手能照着复现，熟手能直接看边界条件。
2. 从波形到标签：地震数据集怎么切、怎么标、怎么存
2.1 地震波形数据的来源与筛选标准
做地震等级预测，第一步不是搭网络，是把数据搞干净。常见来源有三类：IRIS 等公开台网下载的 miniSEED，本地台网归档的连续波形，以及已经切好的事件波形库。我一般优先用事件波形库，因为连续波形切事件这一步本身就容易引入偏差——STA/LTA 触发阈值设高了漏掉小震，设低了把噪声当事件。
筛选标准要卡死几条。第一，震中距控制在 0 到 200 km，太远的高频衰减严重，P 波信噪比撑不住。第二，震级范围按你的目标定，如果做 3.0 到 6.0 的等级预测，两端各留 0.2 的缓冲，避免边界样本把回归拉偏。第三，每个事件至少三个台站有记录，单台站样本对 CNN 来说信息太单薄。第四，采样率统一重采样到 100 Hz，不统一的话卷积核感受野对应的物理时间就不一致，训练出来的模型换个台网就废。
% 读取 miniSEED 并做基本筛选 cfg = struct('fs', 100, 'distMax', 200, 'magMin', 3.0, 'magMax', 6.0); evtList = dir(fullfile(dataRoot, '*.mseed')); validEvt = {}; for k = 1:numel(evtList) [w, meta] = readseis(fullfile(dataRoot, evtList(k).name)); if meta.dist > cfg.distMax || meta.mag < cfg.magMin || meta.mag > cfg.magMax continue; end if meta.fs ~= cfg.fs w = resample(w, cfg.fs, meta.fs); end validEvt{end+1} = struct('wave', w, 'mag', meta.mag, 'dist', meta.dist); end这段逻辑的核心是先把物理约束卡住，再做重采样。distMax和震级范围是硬门槛，resample放在筛选之后是为了省算力。注意readseis是示意函数名，实际用rdsac、irisFetch或自己写的 miniSEED 解析都行，关键是拿到波形矩阵和元数据。
2.2 三分量切窗与标签构造
CNN 的输入张量一般是[时间点 × 通道 × 1]或者[时间点 × 1 × 通道]，MATLAB 深度学习工具箱习惯把通道放最后一维，所以推荐[N × 3 × 1]，N 是窗口长度。窗口怎么切有讲究：以 P 波到时为基准，往前取 1 秒做前置噪声，往后取 9 秒覆盖 S 波和部分尾波，总长 10 秒，100 Hz 下就是 1000 个点。这个 1+9 的分配不是拍脑袋，前置噪声给网络一个信噪比参照，后段尾波携带震级信息。
标签构造有两种路线。一是直接回归震级数值，输出一个标量，损失用 MSE。二是把震级离散成等级，比如 3.0-3.9、4.0-4.9、5.0-5.9 三档，做分类。回归对连续预测更友好，分类在样本不均衡时更稳。我一般先做回归，因为地震等级本质是连续量，硬切等级会丢信息。
winPre = 1 * cfg.fs; % P 波前 1 秒 winPost = 9 * cfg.fs; % P 波后 9 秒 X = zeros(numel(validEvt), winPre+winPost, 3, 'single'); Y = zeros(numel(validEvt), 1, 'single'); for i = 1:numel(validEvt) pIdx = validEvt{i}.pIndex; seg = validEvt{i}.wave(pIdx-winPre+1 : pIdx+winPost, :); seg = seg - mean(seg, 1); % 去均值 seg = seg ./ (max(abs(seg), [], 1) + eps); % 逐通道归一化 X(i,:,:) = seg; Y(i) = validEvt{i}.mag; end去均值和逐通道归一化是必须的，不同台站的仪器响应和增益差异很大，不归一化网络会先学台站特征而不是震源特征。eps防止除零。归一化用max(abs())而不是标准差，是因为地震波形有尖脉冲，标准差会被异常值拉大。
2.3 训练集、验证集、测试集怎么分才不泄漏
随机划分是新手最容易翻车的地方。同一个地震事件被多个台站记录，如果按样本随机分，同一事件的不同台站样本可能同时进训练集和测试集，模型等于见过答案，测试精度虚高。正确做法是按事件 ID 分组划分，同一事件的所有台站样本要么全在训练集，要么全在测试集。
evtIds = unique({validEvt.eventId}); rng(42); idx = randperm(numel(evtIds)); nTrain = round(0.7 * numel(evtIds)); nVal = round(0.15 * numel(evtIds)); trainEvt = evtIds(idx(1:nTrain)); valEvt = evtIds(idx(nTrain+1:nTrain+nVal)); testEvt = evtIds(idx(nTrain+nVal+1:end));比例我一般用 70/15/15。验证集用来早停和调参，测试集只在最后跑一次。rng(42)固定随机种子，保证复现。如果事件数少于 500，建议用 5 折交叉验证代替单次划分，否则测试集太小，指标波动大。
3. 在 MATLAB 里搭一个能用的 CNN 回归网络
3.1 网络结构选型：为什么不用现成的 ResNet
地震波形是一维时序，直接套二维图像网络（ResNet、VGG）是常见误用。二维卷积核在图像上抓的是空间局部性，地震波形只有时间一维，硬把[1000×3]当图像处理，卷积核会在通道维度上乱滑，学到的特征没有物理意义。正确做法是用convolution1dLayer，卷积核只在时间轴上滑动。
网络深度不用太深。我一般用 4 到 6 个卷积块，每块是conv1d + batchNorm + relu + maxPool。第一层卷积核设大一点，比如 15 到 31，感受野覆盖 P 波初动的几个周期；后面逐层减小到 3 到 5，抓更细的尾波衰减模式。通道数从 32 起步，每池化一次翻倍，到 256 封顶。再深就容易过拟合，地震样本量通常撑不住。
layers = [ sequenceInputLayer(3, 'Name', 'input') % 3 通道输入 convolution1dLayer(31, 32, 'Padding', 'same', 'Name', 'conv1') batchNormalizationLayer('Name', 'bn1') reluLayer('Name', 'relu1') maxPooling1dLayer(4, 'Stride', 4, 'Name', 'pool1') convolution1dLayer(15, 64, 'Padding', 'same', 'Name', 'conv2') batchNormalizationLayer('Name', 'bn2') reluLayer('Name', 'relu2') maxPooling1dLayer(4, 'Stride', 4, 'Name', 'pool2') convolution1dLayer(7, 128, 'Padding', 'same', 'Name', 'conv3') batchNormalizationLayer('Name', 'bn3') reluLayer('Name', 'relu3') maxPooling1dLayer(4, 'Stride', 4, 'Name', 'pool3') convolution1dLayer(3, 256, 'Padding', 'same', 'Name', 'conv4') batchNormalizationLayer('Name', 'bn4') reluLayer('Name', 'relu4') globalAveragePooling1dLayer('Name', 'gap') dropoutLayer(0.3, 'Name', 'drop') fullyConnectedLayer(64, 'Name', 'fc1') reluLayer('Name', 'relu5') fullyConnectedLayer(1, 'Name', 'fcOut') % 回归输出标量 regressionLayer('Name', 'reg') ];sequenceInputLayer(3)对应三分量。Padding='same'保证卷积不改变时间长度，池化负责降维。globalAveragePooling1dLayer替代全连接展平，参数量小很多，抗过拟合。最后fullyConnectedLayer(1)加regressionLayer是回归标配。如果做分类，把最后两层换成fullyConnectedLayer(numClasses)和classificationLayer。
3.2 训练参数怎么设：学习率、批大小、早停
训练用trainingOptions配。优化器选adam，学习率初始 1e-3，每 10 个 epoch 降一半。批大小 32 到 64，取决于显存，地震样本 1000 点不算大，64 一般够。epoch 上限设 100，但一定要开验证早停，ValidationPatience设 10，验证损失 10 轮不降就停。
opts = trainingOptions('adam', ... 'InitialLearnRate', 1e-3, ... 'LearnRateSchedule', 'piecewise', ... 'LearnRateDropPeriod', 10, ... 'LearnRateDropFactor', 0.5, ... 'MaxEpochs', 100, ... 'MiniBatchSize', 64, ... 'Shuffle', 'every-epoch', ... 'ValidationData', {XVal, YVal}, ... 'ValidationFrequency', 20, ... 'ValidationPatience', 10, ... 'OutputNetwork', 'best-validation-loss', ... 'Plots', 'training-progress', ... 'Verbose', false);OutputNetwork='best-validation-loss'是关键，训练结束时返回验证损失最低的那版网络，而不是最后一版。Shuffle='every-epoch'打乱样本顺序，避免批次间相关性。学习率衰减用 piecewise，比 cosine 更直观，调起来好控制。
3.3 训练、验证与指标解读
net = trainNetwork(XTrain, YTrain, layers, opts); YPred = predict(net, XTest); rmse = sqrt(mean((YPred - YTest).^2)); mae = mean(abs(YPred - YTest)); r2 = 1 - sum((YTest - YPred).^2) / sum((YTest - mean(YTest)).^2); fprintf('RMSE=%.3f MAE=%.3f R2=%.3f\n', rmse, mae, r2);震级预测里 RMSE 到 0.3 以下算能用，0.2 以下算好。R2 低于 0.7 说明模型没学到东西，先查数据泄漏和归一化。MAE 比 RMSE 抗异常值，两个一起看。如果 RMSE 正常但 MAE 很大，说明有个别样本预测偏得离谱，回去查那些样本的 P 波拾取是不是错了。
4. 避坑与排查：地震 CNN 训练里最容易翻车的五件事
4.1 损失不下降，先查标签和归一化
现象：训练几个 epoch 后 loss 卡在某个值不动，验证 loss 还往上走。原因通常是标签没对齐或者归一化做错了。比如 P 波到时索引偏移了一位，切出来的窗口全是噪声；或者归一化时对整个数据集算了一个全局均值，而不是逐样本逐通道。解决：随机抽 10 个样本画波形图，肉眼确认 P 波在窗口内的位置；归一化改成逐样本逐通道，代码里mean(seg,1)的维度别写错。
4.2 验证集精度远高于测试集
现象：验证 RMSE 0.15，测试 RMSE 0.45。原因几乎都是数据泄漏——按样本随机划分，同一事件的不同台站样本跨了集合。解决：改成按事件 ID 分组划分，代码见 2.3 节。如果已经分组还这样，检查验证集和测试集的事件有没有重叠，用intersect查一下。
4.3 模型对远震预测系统性偏大
现象：震中距大于 150 km 的样本，预测震级普遍比真实值高 0.3 到 0.5。原因是远震高频衰减，波形幅度小，归一化后噪声被放大，网络把噪声当成了大震的尾波。解决：在输入里加一个震中距通道，让网络知道距离信息；或者对远震样本单独做增益补偿。我一般加距离通道，简单有效。
4.4 训练到一半 loss 突然变 NaN
现象：第 30 个 epoch 左右 loss 变 NaN，训练中断。原因是学习率太大，梯度爆炸。解决：把初始学习率降到 1e-4，加梯度裁剪。MATLAB 的trainingOptions没有直接的梯度裁剪参数，可以在自定义训练循环里用dlupdate裁剪，或者干脆把学习率调小、批大小调大。
4.5 换个台网模型就废
现象：在 A 台网训练 R2 0.85，拿到 B 台网测试 R2 掉到 0.4。原因是仪器响应和台基条件不同，模型学到了台站特征。解决：训练时做台站增广，随机给波形乘一个 0.8 到 1.2 的增益因子，模拟不同台站响应；或者用更多台网的数据一起训练，让网络学会忽略台站差异。
5. 把震级预测精度再压 0.05：几个我常用的进阶技巧
第一个技巧是加物理约束特征。纯波形输入网络要自己学 P 波到时，但到时拾取本身有成熟算法，把 STA/LTA 触发点、P 波初动极性、S-P 时差这三个特征拼到全连接层前面，比让卷积核硬学要快。具体做法是在globalAveragePooling1dLayer后面接一个concatenationLayer，把三个标量特征和池化输出拼起来再进全连接。
% 在 gap 后拼接物理特征 featLayer = [ globalAveragePooling1dLayer('Name','gap') concatenationLayer(1, 2, 'Name', 'concat') fullyConnectedLayer(64, 'Name', 'fc1') reluLayer('Name', 'relu5') fullyConnectedLayer(1, 'Name', 'fcOut') regressionLayer('Name', 'reg') ]; % 需要把 layers 里 gap 之后的部分替换掉，并用 dlnetwork 手动前向拼接层要求两个输入维度匹配，物理特征要归一化到和池化输出同量级。这个改动在样本量大于 2000 时提升明显，小于 1000 时可能过拟合，慎用。
第二个技巧是测试时增强（TTA）。对同一个测试样本，做三次微小扰动——时间轴平移 ±50 ms、加 1% 高斯噪声——分别预测，取三次均值。这个操作不增加训练成本，测试阶段多跑两遍，RMSE 一般能降 0.02 到 0.04。代价是推理时间翻三倍，实时场景要权衡。
| 技巧 | 预期 RMSE 降幅 | 适用样本量 | 推理开销 |
|---|---|---|---|
| 物理特征拼接 | 0.03-0.06 | >2000 | 无 |
| 测试时增强 | 0.02-0.04 | 任意 | 3 倍 |
| 多台网联合训练 | 0.05-0.10 | >5000 | 无 |
| 震级分层加权采样 | 0.02-0.05 | 不均衡时 | 无 |
第三个技巧是震级分层加权采样。地震样本天然不均衡，小震多、大震少，回归模型会偏向预测中间值。在trainingOptions里没有直接的样本权重参数，我一般用自定义训练循环，在计算损失时给大震样本乘一个权重，权重取1/sqrt(freq)，freq 是该震级档的样本频率。这个改动对 5.0 以上样本的预测精度提升最明显。
最后一个习惯：每次训练完，把测试集预测值和真实值画散点图，看残差分布。如果残差在某个震级段系统性偏正或偏负，说明模型在那个段有偏差，回去查那段样本的 P 波拾取质量。我踩过最深的坑就是一批 4.5 级样本的 P 波拾取用了错误的到时，模型怎么调都差 0.1，换掉那批标签后直接达标。数据质量永远比网络结构重要，希望帮到你。
本文还有配套的精品资源，点击获取
