简介基于MATLAB编程的卷积神经网络CNN地震等级预测项目面向地震信号分析与深度学习入门者目标是将CNN模型用于地震等级预测解决传统信号处理方法特征提取依赖手工设计的问题。资源共33个文件以30个m脚本为主涵盖网络搭建、参数初始化、前向/反向传播、梯度检查、超参数组合遍历、数据预处理等环节配合2个mat数据/模型文件和1个xlsx表格压缩包仅370KB结构紧凑便于整体研读。已有350人学习适合需要通过完整代码理解CNN底层计算过程的读者。通过这些代码读者能复现从地震数据读入、特征变换到模型训练与预测的完整流程也能通过逐行阅读脚本掌握卷积层、池化层、全连接层及损失计算的具体实现为后续调整网络结构或迁移到其他信号分类任务提供可运行的起点。1. 地震等级预测为什么选 CNN 而不是传统回归一次中等地震发生后不同台站基于波形振幅和周期估算的震级常常相差 0.3 级以上值班人员要依靠多台平均、人工剔异才能给出稳定数值。这个过程的本质是从波形到标量的回归给出一段三分量记录让它对应一个连续的震级值M_L 或 M_w而不是让模型去“识别”地震大小。基于MATALB编程的卷积神经网络CNN处理这件事是把波形当作一维时间序列让卷积层自动提取震相触发、振幅包络和衰减特征端到端输出震级。相比人工挑特征的传统回归管线CNN 的优势在样本量上来之后才明显适合台网密集、事件目录完整但要反复为每个新区块做场地校正的场景。下面按数据准备、网络搭建、训练调参、验证修正四个部分展开每个环节都给出能直接跑的最小方案。2. 地震数据预处理波形对齐、归一化与标签构造2.1 数据从哪来怎么组织成训练样本训练 CNN 需要的是成对数据台站记录到的波形片段配上该次地震的震级标签。数据来源一般是公开地震波形库如 IRIS、GFZ 提供的 SAC 或 miniseed 格式数据也可以是自己台网中心导出的归档波形。关键不在格式而在每条记录必须带上三个元信息事件编号、台站编号、P 波或 S 波到时。没有到时标注后续截窗就是空的。样本组织建议采用“一个台站对一次地震的一条记录”为一条样本。三分量可以先合并成单通道如合成矢量幅值也可以保留三个分量作为三通道输入。前者实现简单后者信息更全但输入维度增加后网络参数和训练时间都会上升。给出一张组织格式表便于对照你手头的数据字段含义示例波形数组某台站某次地震的采样序列4000 个点 40 Hz事件编号对应地震目录里的唯一 ID20230042台站编号用于分组和后期残差修正ST01标签该地震的目录震级4.6P 波到时用于截窗对齐第 1050 个采样点这里有一个容易忽略的点如果你要估计的是面波震级 M_s 或矩震级 M_w波形窗口必须包含尾波和面波段只截 P 波到后几秒是远远不够的。常见做法是截取 P 波到时前 10 秒、到时后 90 秒这样模型能看到完整的震相序列和衰减过程。2.2 预处理步骤与最小可运行代码预处理通常包含四步去均值、幅值归一化、重采样、按 P 波到时截窗。去均值是去掉记录里的直流漂移幅值归一化是为了消除台站增益差异带来的量纲问题重采样是为了统一采样率并降低序列长度截窗则是把波形送到网络之前的最后一步对齐。下面这段是 MATALB 中等价于“标准四步”的最小实现% data: N x T 原始波形矩阵已按事件-台站展开 % Fs0: 原始采样率单位 Hz % labels: N x 1 震级标签 Fs0 100; % 假设原始采样率 100 Hz fs 40; % 目标采样率 40 Hz保留 0.1~15 Hz 主频带 pre 10; % P 波到时前保留 10 s post 90; % P 波到时后保留 90 s L (pre post) * fs; X zeros(size(data,1), L); for i 1:size(data,1) x data(i,:); x x - mean(x); % 去直流 x x / (max(abs(x)) eps); % 幅值归一化防止台站增益差异主导 xr resample(double(x), fs, Fs0); % 重采样内部自带抗混叠滤波 % 假设每条记录已经把 P 波到时放在序列中点 c0 round((pre post)/2 * fs); seg xr(c0 - pre*fs : c0 post*fs - 1); seg seg / (std(seg) eps); % 再按能量归一化一次 X(i,:) seg; end这段代码有两个参数值得说明。pre10和post90决定了网络能看到 P 波前 10 秒的环境噪声和 P 波后 90 秒的完整衰减如果目标震级类型换成近震震级 M_L窗口可以缩到 30 秒以内训练速度会明显更快。resample来自 Signal Processing Toolbox函数内部会做抗混叠滤波所以不需要在它前面再串一个巴特沃斯低通。第二次归一化用std而不是最大值是为了防止单个尖锐脉冲把整段能量压得过低这在有脉冲干扰的记录里很常见。2.3 数据泄漏按事件划分训练集和测试集这是新手最容易踩的一个坑也是检验预处理环节是否专业的分水岭划分训练集和测试集时必须按事件划分而不是按记录划分。同一个地震事件往往有几十个台站同时记录到如果随机打乱记录同一次地震的台站记录会同时出现在训练集和测试集中模型相当于在考场上见过标准答案验证指标会虚高 20% 甚至更多。按事件划分的做法是先取出全部事件编号随机抽 20% 的事件作为测试事件再把这些事件对应的所有台站记录划入测试集rng(42); allEvents unique(eventId); nTest max(1, round(numel(allEvents) * 0.2)); testEvents allEvents(randperm(numel(allEvents), nTest)); isTest ismember(eventId, testEvents); XTrain X(~isTest, :); YTrain labels(~isTest, :); XTest X(isTest, :); YTest labels(isTest, :);注意这里的eventId是每条记录对应的事件编号不是台站编号。划分完成后可以检查一下isTest中是否包含完整的事件组避免因为编号排序问题导致同一个事件的记录被切成两半。这样划分出来的测试集效果才接近真实台网部署场景一条新记录到来时模型之前从没见过这个事件。3. 用 MATALB 搭建 CNN卷积层参数与网络结构3.1 输入形态决定用 1D-CNN 还是 2D-CNN很多第一次做地震波形的同学会把波形转成语谱图去套图像分类网络这其实是把简单问题复杂化。转换成二维语谱图需要额外决定窗长、窗移、频率分辨率等一系列超参数等于把特征工程的负担从网络结构转移到了预处理阶段而且引入的时频表示会丢失原始波形的相位信息。地震波形本质上是时间序列直接用一维卷积作用于原始波形是最直接的做法模型需要什么频率特征卷积核自己会学。在 MATALB 里对应的是convolution1dLayer输入层用sequenceInputLayer。如果只用单通道的合成幅度波形第一层写sequenceInputLayer(1)如果想把垂直分量、东西分量、南北分量作为三通道写sequenceInputLayer(3)。多通道输入不会改变整体结构只是让第一次卷积同时看到三个分量的信息在样本量够大时通常有收益。3.2 最小化 CNN 结构代码与参数说明下面是一个可以直接复制到 Deep Learning Toolbox 里跑的回归型 CNN 结构。它由两个卷积块加全局平均池化组成输出层只有一个神经元对应连续的震级值。regressionLayer会自动使用均方误差作为损失函数。layers [ sequenceInputLayer(1, Name, input) convolution1dLayer(20, 16, Padding, same, Name, conv1) batchNormalizationLayer(Name, bn1) reluLayer(Name, relu1) maxPooling1dLayer(4, Stride, 4, Name, pool1) convolution1dLayer(12, 32, Padding, same, Name, conv2) batchNormalizationLayer(Name, bn2) reluLayer(Name, relu2) maxPooling1dLayer(4, Stride, 4, Name, pool2) globalAveragePooling1dLayer(Name, gap) fullyConnectedLayer(1, Name, fc) regressionLayer(Name, output) ]; lgraph layerGraph(layers); analyzeNetwork(lgraph); % 可生成 CNN 结构图和每一层的尺寸信息用analyzeNetwork弹出来的窗口可以检查每一层的输出维度这是排查输入尺寸不匹配最直接的工具。参数层面convolution1dLayer(20, 16)的第一个参数 20 是卷积核长度40 Hz 采样率下对应 0.5 秒的时间跨度刚好覆盖一个完整的震相周期第二个参数 16 是输出通道数。maxPooling1dLayer(4, Stride, 4)表示每 4 个点取一次最大值效果是把时间轴分辨率降低四倍让后续层看到更大的时间范围。3.3 感受野计算与网络深度卷积网络的一个核心问题是感受野最后一层特征上的一个点对应原始波形上的多少个采样点。这个数值决定了模型能否“看到”完整的震相序列。感受野的递推计算方式是从输入层往后每经过一层当前感受野乘以该层的步长再加上该层卷积核带来的增量。以刚才的结构为例层核长/步长感受野增量累计感受野conv120 / 120 点20 点pool14 / 4020 点conv212 / 1(12-1) × 4 44 点64 点pool24 / 4064 点64 点在 40 Hz 采样率下只有 1.6 秒。也就是说上面这个网络实际只能看到 P 波到达后约 1.6 秒内的波形对震级估计来说远远不够。要覆盖 15 秒以上的波形常见做法是把网络加深到三个或四个卷积块同时保持池化层步长为 4 快速下采样。加深之后的感受野大约可以到几百点但这会带来参数量的膨胀和训练速度的下降。如果不想牺牲结构可以在第 5 章用多尺度输入的方式绕开感受野限制这个后面会具体说。4. 训练 CNN 地震预测模型损失、学习率与验证策略4.1 回归与分类损失函数怎么选标题里的“地震等级”在某些语境下被当成分类问题处理这需要先澄清。震级本身是连续值M_L 4.6 和 M_L 4.7 之间的差别是有物理意义的用回归损失更合理。如果模型输出层接的是classificationLayer训练目标变成“这个震级属于 4.0~4.9 还是 5.0~5.9”精度会受分箱边界影响而且震级落在边界附近时分类结果极不稳定。所以这里选regressionLayer损失函数为均方误差。如果确实要做烈度预测这种序数分类任务不建议直接套分类网络而是把输出层改为一个神经元仍然用回归把烈度等级编码为数值预测后取最近整数。这样保留了等级之间的顺序关系训练也更稳定。4.2 训练选项的关键参数训练回归网络使用的优化器通常是adam它对初始学习率不敏感是地震这类高噪声数据最不容易发散的默认选择。下面是训练代码和一份带解释的参数表options trainingOptions(adam, ... MaxEpochs, 60, ... MiniBatchSize, 32, ... InitialLearnRate, 2e-3, ... LearnRateSchedule, piecewise, ... LearnRateDropFactor, 0.3, ... LearnRateDropPeriod, 20, ... L2Regularization, 5e-4, ... ValidationData, {XVal, YVal}, ... ValidationFrequency, 10, ... Shuffle, every-epoch, ... Plots, training-progress); net trainNetwork(XTrain, YTrain, lgraph, options);参数建议值说明InitialLearnRate1e-3 ~ 3e-3大于 1e-2 时 loss 很容易发散成 NaNMiniBatchSize16 ~ 64由显存决定OOM 时优先减半LearnRateDropFactor0.3每 20 轮降低为原来的 30%后期微调权重ValidationFrequency10每 10 轮在验证集上评估一次L2Regularization5e-4卷积层参数多正则系数过大容易欠拟合这里有一个训练技巧值得单独提InitialLearnRate不能只看 loss 曲线的起点还要观察前 5 个 epoch 内验证集的表现。如果验证 loss 在第 3 轮左右突然跳高多半是学习率偏大如果整个训练过程验证 loss 都不下降则可能是这个学习率太小需要从 1e-2 起做几次 warmup 找到合适区间。现在很多同事习惯让 AI 编程辅助生成训练骨架但训练参数这部分必须自己逐项确认AI 生成的超参组合通常来自图像分类实验直接用到波形回归上并不合适。4.3 常见训练故障与排查训练地震波形 CNN 最容易遇到三类问题。第一类是 loss 直接变 NaN原因通常是输入数据里存在 NaN 值或者学习率过大先检查 X 里有没有isnan再调低学习率。第二类是训练集 loss 持续下降但验证集 loss 在某个 epoch 后反弹这是典型的过拟合信号优先减小网络通道数或提高 L2 正则。第三类问题比较隐蔽验证集 loss 看起来很低但把真实震级和预测值画成散点后斜率明显小于 1说明模型学会了输出平均值附近的值这是回归问题中 MSE 损失常见的“回归到均值”现象。应对方式是检查数据集中震级的分布是否过于集中如果大部分样本在 4.0~5.0 之间考虑按震级分层采样或者改用 Huber 损失来提高对极端样本的敏感度。5. 验证技巧多尺度输入与台站残差修正5.1 多尺度时间窗让 CNN 同时看到 P 波和尾波前面说过加深网络可以扩大感受野但代价是网络更重、训练更慢。另一个更轻量的替代方案是多尺度输入把同一个事件截取成不同窗口长度分别重采样到相同点数再叠成多通道输入。短窗口保留 P 波的细节触发特征长窗口提供尾波衰减的整体信息。实现时只需要改输入层通道数网络结构不变% 构造三通道输入15s / 30s / 90s 三种窗口 ch zeros(3, 4000); ch(1,:) extractWindow(x, 15, fs); ch(2,:) extractWindow(x, 30, fs); ch(3,:) extractWindow(x, 90, fs); % 网络第一层改为 sequenceInputLayer(3, Name, input3)需要说明的是三个窗口长度对应的时间范围不同但长度经过重采样后都统一成 4000 点这样模型会把每个通道当作不同的“视角”而不是不同的长度。训练时对每个通道分别做归一化避免短窗口的高频细节被长窗口的能量淹没。5.2 台站残差修正与震级一致性检查如果说多尺度输入解决的是信号信息量问题台站残差修正解决的是系统偏差问题。每个台站由于台基条件、仪器响应和场地效应的差异对同一次地震的震级估计会有固定的系统偏差这个偏差不会因为网络结构变化而消失。做法是在训练集上预测一遍按台站编号计算预测值与真实值的平均残差之后对该台站的新预测做修正resid predict(net, XTrain) - YTrain; bias accumarray(stationId, resid, [], mean); predCorrected predRaw - bias(stationIdNew);这个修正操作相当于在模型输出上叠加了一层台站级校准效果比调整网络结构更直接。实践中我一般会在台站偏差计算之后画一张偏差随方位角变化的散点图如果偏差呈现明显的随方位角振荡说明还有三维介质各向异性在起作用这时仅靠 CNN 无法完全建模可以考虑把震中距作为额外特征拼接进全连接层。这样修完的预测值才能进入后续的目录产出环节。本文还有配套的精品资源点击获取
