简介:基于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 |
| 事件编号 | 对应地震目录里的唯一 ID | 20230042 |
| 台站编号 | 用于分组和后期残差修正 | ST01 |
| 标签 | 该地震的目录震级 | 4.6 |
| P 波到时 | 用于截窗对齐 | 第 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这段代码有两个参数值得说明。pre=10和post=90决定了网络能看到 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 感受野计算与网络深度
卷积网络的一个核心问题是感受野:最后一层特征上的一个点,对应原始波形上的多少个采样点。这个数值决定了模型能否“看到”完整的震相序列。感受野的递推计算方式是,从输入层往后,每经过一层,当前感受野乘以该层的步长,再加上该层卷积核带来的增量。以刚才的结构为例:
| 层 | 核长/步长 | 感受野增量 | 累计感受野 |
|---|---|---|---|
| conv1 | 20 / 1 | 20 点 | 20 点 |
| pool1 | 4 / 4 | 0 | 20 点 |
| conv2 | 12 / 1 | (12-1) × 4 = 44 点 | 64 点 |
| pool2 | 4 / 4 | 0 | 64 点 |
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);| 参数 | 建议值 | 说明 |
|---|---|---|
| InitialLearnRate | 1e-3 ~ 3e-3 | 大于 1e-2 时 loss 很容易发散成 NaN |
| MiniBatchSize | 16 ~ 64 | 由显存决定,OOM 时优先减半 |
| LearnRateDropFactor | 0.3 | 每 20 轮降低为原来的 30%,后期微调权重 |
| ValidationFrequency | 10 | 每 10 轮在验证集上评估一次 |
| L2Regularization | 5e-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 无法完全建模,可以考虑把震中距作为额外特征拼接进全连接层。这样修完的预测值,才能进入后续的目录产出环节。
本文还有配套的精品资源,点击获取