简介：使用LSTM神经网络进行地震震级预测，是地震数据分析中典型的深度学习应用场景。本资源即围绕这一主题，提供一套基于Matlab实现的完整工程代码，适用于计算机、电子信息工程、数学等专业学生在课程设计、期末大作业或毕业设计中快速上手。整个压缩包共49个文件，容量仅1.35MB，包含Matlab脚本（.m）、Python辅助脚本（.py）、CSV数据文件、PNG结果图及FIG图形窗口等，既保留了核心预测算法，也附带了丰富的可视化分析结果，便于对照学习。目前已有176人学习下载。代码采用参数化编程，思路清晰、注释详细，并内置可直接运行的案例数据，可帮助读者快速复现地震震级预测流程。资源中不仅包含LSTM时序预测模型，还提供了数据预处理、可视化工具以及多张分析图表，如震级与深度散点图、震级直方图、地震数据时间序列等，能够支撑从数据读取、模型训练到结果评估的完整链路，是理解LSTM在地震领域应用的实用参考资料。

1. 为什么地震序列适合 LSTM：从震级预测到波形分析

地震预测是个老话题，但地震震级预测直到今天仍是个不确定性问题。传统做法多依赖统计模型和物理模型，比如 Gutenberg-Richter 定律做长期概率估计，或者用库仑应力变化做短期触发判断。这些方法的问题在于，它们很难利用地震目录里蕴含的时间顺序信息——余震序列的衰减规律、主震前的小震群活动模式、不同震级区间的演化节奏，这些特征本质上都是时间序列特征，适合用循环神经网络来建模。LSTM 作为循环神经网络的一种变体，通过门控机制解决长序列训练中的梯度消失问题，能够在一定时间窗口内记住前序地震事件的震级、时间和空间信息，从而对下一次可能发生的震级做出估计。

本标题里的「matlab代码.zip」表明这套方案不需要额外安装 PyTorch 或 TensorFlow，直接使用 MATLAB 的 Deep Learning Toolbox 就能完成数据读取、网络训练和震级预测全过程。适合三类人：一是有地震目录数据但缺乏深度学习落地经验的研究生和工程师；二是想用 MATLAB 完成时序预测但对 LSTM 参数调节不熟的地球物理从业者；三是想快速验证 LSTM 在震级预测上效果、又不想陷入 Python 环境配置的算法工程师。接下来的内容会从数据预处理、模型构建、参数设置到评估验证逐步展开，中间会给出能直接改路径运行的 MATLAB 代码。

2. 地震目录数据的预处理：从 Excel 台账到适合 LSTM 输入的序列样本

2.1 理解 LSTM 的输入格式要求

LSTM 网络在 MATLAB 中的输入通常是一个元胞数组形式的序列，每个元素对应一个时间步的输入特征向量。对于地震震级预测来说，特征向量可以是单一震级标量，但更有效的做法是把时间间隔、震级、纬度、经度、深度以及历史最大震级一起作为特征。原因在于，地震目录中每次事件都有区域信息，如果只把震级序列喂给 LSTM，模型实际上只能学到震级本身的统计规律，而无法捕捉同一断层带或同一网格内地震的空间依赖关系。

常见做法是将目录数据按时间排序后，做分箱处理。比如按日或按小时把研究区域划分成 0.1°×0.1° 的网格，每个网格作为一条独立序列，序列中每个时间点的特征向量包含该网格内是否发生地震、最大震级、平均震级、累计能量对数等。另一种更简单但有效的做法是只取全局地震目录中最大震级高于阈值的事件，构建一条全局震级序列。对于震级预测这个任务，我一般先做全局序列，因为模型复杂度低、调试快，效果好再加空间特征。

% 读取地震目录（假设 CSV 包含 Date, Time, Lat, Lon, Depth, Magnitude 列）
data = readtable('earthquake_catalog.csv');
t = datetime(data.Date, 'InputFormat', 'yyyy-MM-dd');
t = t + duration(data.Time); % 合并日期和时间

% 按时间排序
[~, idx] = sort(t);
data = data(idx, :);
t = t(idx);

% 筛选最小震级，避免背景小震干扰
minMag = 3.0;
data = data(data.Magnitude >= minMag, :);
t = t(data.Magnitude >= minMag);

disp(['有效事件数: ', num2str(height(data))]);

这段代码的核心是把原始目录整理成一条按时间排列的事件序列。 readtable 读取的 CSV 需要包含至少五列，日期时间解析时要注意 MATLAB 的 datetime 格式字符串必须与 CSV 中日期列的实际格式严格匹配，否则会出现 NaT 值。筛选最小震级的目的是减少噪声事件，因为小震级事件大量存在，会让 LSTM 把精力分散在学习背景噪声上，而不是主要震级变化趋势。

2.2 构造监督学习样本：用时间窗口切出训练集

LSTM 预测震级通常采用滑窗方式构造样本。假设使用前 10 个连续事件窗口来预测下一个事件震级，每个样本的输入是 10×F 的矩阵（F 为特征维度），输出是第 11 个事件的震级。滑窗之后，样本之间的重叠会造成一定程度的数据泄漏，但在时序预测中这是标准做法，关键是保证训练集和测试集按时间切分，不能让测试集的样本混入训练集。

windowSize = 10;      % 用过去10个事件预测下一个
numFeatures = 4;      % 特征数：震级、时间差、深度、纬度
X = {};
Y = [];

for i = 1:(height(data) - windowSize)
    % 构建时间差：当前事件与上一个事件的小时间隔小时
    dt = hours(t(i+1:windowSize+i) - t(i:windowSize+i-1));
    % 特征矩阵
    seq = [data.Magnitude(i:i+windowSize-1), ...
           dt, ...
           data.Depth(i:i+windowSize-1), ...
           data.Lat(i:i+windowSize-1)];
    X{end+1, 1} = seq';  % LSTM要求序列维度：特征数×窗口长度
    Y(end+1, 1) = data.Magnitude(i+windowSize);
end

% 按时间顺序划分训练/验证/测试集（不随机打乱）
numTrain = floor(0.7 * length(X));
numVal = floor(0.15 * length(X));
XTrain = X(1:numTrain);
YTrain = Y(1:numTrain);
XVal = X(numTrain+1:numTrain+numVal);
YVal = Y(numTrain+1:numTrain+numVal);
XTest = X(numTrain+numVal+1:end);
YTest = Y(numTrain+numVal+1:end);

这里的 seq' 转置很关键：MATLAB 的 LSTM 层要求输入是 numFeatures × sequenceLength 的矩阵，与 Python 中 (batch, seq_len, features) 的格式完全不同，初次从 PyTorch 转过来的人最容易在这里报错。 dt 时间差用小时为单位，目的是让 LSTM 感知到事件之间的时间间隔——余震序列在初期密集、后期稀疏，时间差本身就是一个强特征，对于判断当前处于主震后的哪个阶段非常有价值。深度和纬度同样参与计算，相当于给每个事件一个空间坐标上下文。

2.3 归一化：震级特征必须做，但要注意反归一化的位置

LSTM 内部使用 tanh 和 sigmoid 激活函数，输入特征如果范围差异太大（比如纬度是 30~40，而震级是 3~7），会导致梯度更新过程中某些门控信号饱和，训练不稳定。因此输入特征必须归一化到接近 [-1, 1] 或 [0, 1] 的范围。

% 合并所有特征做 min-max 归一化
allData = cell2mat(X);
xmin = min(allData, [], 2);
xmax = max(allData, [], 2);

XNormalized = cellfun(@(seq) (seq - xmin) ./ (xmax - xmin), X, 'UniformOutput', false);
XTrain = XNormalized(1:numTrain);
XVal = XNormalized(numTrain+1:numTrain+numVal);
XTest = XNormalized(numTrain+numVal+1:end);

% 输出 Y 单独归一化
ymin = min(Y);
ymax = max(Y);
YTrainNorm = (YTrain - ymin) / (ymax - ymin);
YValNorm = (YVal - ymin) / (ymax - ymin);
YTestNorm = (YTest - ymin) / (ymax - ymin);

注意输出 Y 的归一化参数 ymin 和 ymax 只用训练集数据计算，测试集的归一化也必须使用训练集的参数。否则，如果测试集里出现比训练集更大的震级，就会导致归一化后值大于 1，预测时反归一化会得到异常结果。另外， cellfun 中每个序列的归一化是使用全局极值完成的，也就是说不管哪个时间窗，相同的原始震级值归一化后是相同的。这样做的好处是 LSTM 看到的是同一个特征空间中的一致映射，不会因样本不同而变化。

3. 用 MATLAB 构建 LSTM 震级预测网络：层结构、dropout 与训练选项

3.1 网络层的选择：sequenceInputLayer 到 fullyConnectedLayer 的合理搭配

LSTM 网络的层设计没有一个绝对最优解，但针对震级预测这种小样本、单变量或多变量回归任务，我通常会采用两层 LSTM 加一层全连接输出的结构。第一层 LSTM 使用较多的隐藏单元来学习短期依赖（比如余震序列的密集程度），第二层 LSTM 隐藏单元减半，进一步提炼长期依赖模式。这种结构在 MATLAB 中可以很直观地定义。

numHiddenUnits1 = 128;
numHiddenUnits2 = 64;
outputMode = 'last';  % 只取最后一个时间步的输出

layers = [
    sequenceInputLayer(numFeatures)
    
    lstmLayer(numHiddenUnits1, 'OutputMode', 'last')
    dropoutLayer(0.2)
    
    lstmLayer(numHiddenUnits2, 'OutputMode', 'last')
    dropoutLayer(0.2)
    
    fullyConnectedLayer(32)
    reluLayer
    
    fullyConnectedLayer(1)
    regressionLayer
];

OutputMode 参数必须设置正确。因为这里是用 10 个时间步预测第 11 个震级值，整个窗口结束后只需要一个输出，所以第一层和第二层 LSTM 都使用 'last' 模式。如果使用 'sequence' 模式，LSTM 会输出每个时间步的隐状态，后面的 fullyConnectedLayer 就会收到一个序列输出，需要额外处理，且对于震级预测这种单点输出任务是浪费计算量。 dropoutLayer(0.2) 放在每个 LSTM 层后面可以有效减少过拟合，因为地震目录数据量通常不大（数千到数万条），而 LSTM 参数量动辄上万，没有 dropout 很容易记住训练集噪声。

3.2 trainingOptions 的参数范围和选择逻辑

MATLAB 的 trainingOptions 控制训练的全部行为，参数设置有明显的经验法则。对于震级预测，我一般这么设置：

参数名	推荐值	原因
MaxEpochs	200	小数据量需要足够轮次收敛，早期停止会自动截断
MiniBatchSize	64	震级序列样本数量级通常是几千到几万，64 是稳定性和速度的折中
InitialLearnRate	0.005	偏高一点，因为 LSTM 对学习率不敏感，配合 Adam 可较快收敛
ValidationFrequency	20	每 20 个 iteration 验证一次，既能观察趋势又不影响速度
GradientThreshold	1	防止个别样本产生的梯度尖峰炸掉训练
LearnRateSchedule	'piecewise'	训练中段降学习率，细化收敛

options = trainingOptions('adam', ...
    'MaxEpochs', 200, ...
    'MiniBatchSize', 64, ...
    'InitialLearnRate', 0.005, ...
    'ValidationData', {XVal, YValNorm}, ...
    'ValidationFrequency', 30, ...
    'Shuffle', 'every-epoch', ...
    'GradientThreshold', 1, ...
    'LearnRateSchedule', 'piecewise', ...
    'LearnRateDropPeriod', 100, ...
    'LearnRateDropFactor', 0.5, ...
    'Plots', 'training-progress', ...
    'Verbose', true);

Shuffle 这里用的是 'every-epoch' ，它打乱的是每个 epoch 内部 mini-batch 的选取顺序，不会破坏序列内部的时间结构，因为每个样本本身已经是一个完整的时间窗序列。 ValidationData 使用验证集可以让训练过程实时显示验证损失，如果验证损失连续多个 epoch 不降反升，说明过拟合已经开始，这时候应该停止训练，而不是等 200 轮跑完。 GradientThreshold 设为 1 意味着梯度超过 1 就按 1 缩放，这个值在时序预测中几乎可以默认设置，不必轻易改动。

3.3 训练与保存模型：从 trainNetwork 到 matfile 持久化

net = trainNetwork(XTrain, YTrainNorm, layers, options);

% 验证集预测
YPredValNorm = predict(net, XVal, 'MiniBatchSize', 64);
YPredVal = YPredValNorm * (ymax - ymin) + ymin;   % 反归一化

% 保存完整模型
save('lstm_magnitude_model.mat', 'net', 'xmin', 'xmax', 'ymin', 'ymax');

训练完成后， predict 返回的是归一化空间的预测值，必须用之前保存的 ymin 和 ymax 做反归一化。这一步如果漏掉，预测值会全部落在 [0,1] 区间内，看起来完全失真。保存模型时把归一化参数一起存进 .mat 文件是必须的，因为未来加载模型做预测时，新输入数据也要用同样的 xmin / xmax 做归一化处理，才能得到正确的震级输出。

4. 评估震级预测效果：误差指标、波形对比和过拟合判断

4.1 回归任务的三个核心指标

震级预测是个回归问题，常用三个指标来评价：RMSE（均方根误差）、MAE（平均绝对误差）和 R²（决定系数）。RMSE 对大误差更敏感，如果一个预测值与真实值差 1 级，RMSE 会放大这个差异，适合检测极端错误；MAE 则是所有误差的平均水平；R² 说明模型解释了真实震级方差的百分比。在震级预测中，R² 能达到 0.3~0.5 已经算是不错的结果，因为地震本身有很强的随机性，不像股票价格那样有清晰的趋势。

rmse = sqrt(mean((YPredVal - YVal).^2));
mae = mean(abs(YPredVal - YVal));
ssRes = sum((YVal - YPredVal).^2);
ssTot = sum((YVal - mean(YVal)).^2);
r2 = 1 - ssRes / ssTot;

fprintf('RMSE = %.3f 级\n', rmse);
fprintf('MAE  = %.3f 级\n', mae);
fprintf('R²   = %.3f\n', r2);

对于震级预测，RMSE 在 0.5 以下说明模型的预测精度在一个可以讨论的范围内；如果是 0.2 左右，说明模型已经充分学到了训练集中震级变化的规律。但要注意，如果测试集本身震级范围很小（比如测试期间全是 3~3.5 级的小震），RMSE 会天然偏低，这时候要结合 R² 来看——R² 低说明虽然误差绝对值小，但模型并没有解释震级的相对变化，只是预测了平均水平。

4.2 预测曲线对比和错误分布

除了数值指标，可视化是更直观的评估手段。绘制测试集的真实震级与预测震级折线图，可以肉眼看出模型是否跟得上震级突变。LSTM 模型的典型问题是对陡增震级反应迟缓——比如前 5 个事件都是 3.2 级左右，突然发生一个 5.8 级事件，LSTM 往往预测值只在 3.5~4.0 之间。这不是 bug，而是 LSTM 的平滑特性在起作用，它天然会偏向训练集中最常见的震级区间。

figure;
plot(YTest, 'b-', 'LineWidth', 1.2); hold on;
plot(YPredTest, 'r--', 'LineWidth', 1.0);
legend('真实震级', '预测震级');
xlabel('测试集事件序号');
ylabel('震级');
title('LSTM 震级预测 vs 真实震级');
grid on;

错误分布直方图也是一个被忽视但很重要的诊断手段。计算每个样本的预测误差（预测值 - 真实值），如果直方图明显偏斜，模型存在系统性偏差，比如预测值普遍偏高，说明训练集中大震级事件偏多，模型学会了「往大了猜」。此时需要对训练集做震级分桶重采样，让不同震级区间的事件数量大致均匀。

4.3 用验证集判断过拟合和欠拟合

训练过程中的 loss 曲线能反映出训练状态。如果在训练结束时，训练 RMSE 显著低于验证 RMSE，说明过拟合了，应该增大 dropout、减小隐藏单元数或引入 L2 正则化；如果两个 RMSE 都很高且接近，说明模型容量不够，需要增加 LSTM 层数或隐藏单元数。还有一个常见做法是记录每个 epoch 的验证损失，取验证损失最小的那次作为最终模型。在 MATLAB 中， trainNetwork 会自动保存验证损失最小的模型到 net 中，不需要手动判断。

% 预测测试集
YPredTestNorm = predict(net, XTest, 'MiniBatchSize', 64);
YPredTest = YPredTestNorm * (ymax - ymin) + ymin;

% 超过 1 级误差的样本数量
largeErrIdx = abs(YPredTest - YTest) > 1.0;
fprintf('误差超过 1 级的样本占比: %.1f%%\n', 100*mean(largeErrIdx));

这里的 largeErrIdx 统计极端误差样本。如果该比例超过 20%，说明模型在震级骤变场景下表现较差，可以考虑在训练时把相邻事件震级差值作为额外特征输入，或者使用自定义损失函数对大误差施加更大惩罚。

5. 参数调节的三个技巧：窗口长度、隐藏单元数和特征组合的影响

5.1 窗口长度的动态调节策略

windowSize 是最影响预测效果的超参数之一，它决定 LSTM 看到多长的历史。对于地震目录，窗口太短（比如 3~5）会让模型缺乏上下文，无法区分主震后早期余震和背景地震活动；窗口太长（比如 50 以上）则可能引入过多噪声，因为长时间跨度中地震活动的相关性会显著衰减。常见做法是用验证集做一次网格搜索，测试 [5, 10, 15, 20, 30] 几组窗口长度对验证 RMSE 的影响，同时保持其他参数固定。

windowCandidates = [5, 10, 15, 20, 30];
for w = windowCandidates
    % 重新构造样本
    % ...（重复第 2 节的样本构造逻辑）
    % 训练、预测、计算验证 RMSE
    fprintf('Windows = %d, Val RMSE = %.3f\n', w, rmseVal);
end

如果不同窗口长度的效果差异不大，选择最小的值即可，因为更短的窗口意味着样本数更多，模型能获得更多训练数据，并且计算速度更快。如果窗口增大后验证误差明显下降，说明震级序列存在较长的自相关，可以沿用更大窗口，但要注意训练样本数量会成比例减少。

5.2 隐藏单元数的经验边界

隐藏单元数决定了 LSTM 的表示能力，但并不是越大越好。对于小规模地震目录（1000~5000 个事件），我建议隐藏单元数从 32 起步，最高不要超过 256。原因在于，LSTM 的参数量约为 4 * (hiddenSize^2 + hiddenSize*numFeatures) ，128 个隐藏单元意味着第一层 LSTM 就有大约 7 万参数。如果没有足够的数据约束，这些参数会迅速陷入过拟合。手动调节时，把学习率同步调低——隐藏单元翻倍，学习率可以除以 2 或保持不动，防止梯度爆炸。

另一个更有效的做法是单层替代双层。如果数据量不大（小于 3000 个样本），单层 LSTM 加 64 个隐藏单元可能比双层结构更好。因为双层 LSTM 的梯度路径长，小数据量下第二层学不到有意义的隐状态表示。

5.3 特征组合：震级差和时间差的增益

在基础特征（震级、时间差、深度、纬度）之外，我通常会额外构造两个衍生特征：震级差和累积能量。震级差是当前事件震级减上一个事件震级，直接刻画了震级变化的幅度；累积能量用 Gutenberg-Richter 关系的能量公式计算，定义 E = 10^(1.5M + 4.8) ，再把窗口内的能量求和取对数。这个特征可以捕捉到窗口内整体释放能量的趋势，对后续震级预测有额外的信号价值。

% 构造震级差特征
magDiff = [0; diff(data.Magnitude)];

% 构造累计能量特征（窗口内）
energy = 10.^(1.5 * data.Magnitude + 4.8);
logEnergy = zeros(height(data), 1);
for i = 1:height(data)
    startIdx = max(1, i-windowSize+1);
    logEnergy(i) = log10(sum(energy(startIdx:i)));
end

添加特征后， numFeatures 从 4 变为 6，需要同步更新 sequenceInputLayer 的输入维度。这是个常见错误——改动特征矩阵后忘记改网络输入层，训练直接报错。另外，新特征加入后要重新计算 xmin 和 xmax 做归一化，旧归一化参数不能继续使用。

6. 模型部署与持续预测：用训练好的网络实时估算下一个震级

function predMag = predictNextMagnitude(modelFile, recentEvents)
% 加载模型及归一化参数
S = load(modelFile);
net = S.net;
xmin = S.xmin; xmax = S.xmax;
ymin = S.ymin; ymax = S.ymax;

% recentEvents: table，包含 Magnitude, Time, Depth, Lat 列
% 构造特征矩阵，与训练时保持一致
% 注意：时间差需要计算与上一事件的时间间隔
dt = hours(recentEvents.Time(2:end) - recentEvents.Time(1:end-1));
dt = [0; dt];
seq = [recentEvents.Magnitude, dt, recentEvents.Depth, recentEvents.Lat];
seq = seq';

% 归一化
seqNorm = (seq - xmin) ./ (xmax - xmin);

% 预测
predNorm = predict(net, {seqNorm}, 'MiniBatchSize', 1);
predMag = predNorm * (ymax - ymin) + ymin;
end

这个函数的部署逻辑很简单：加载模型和归一化参数，传入最新的一段事件序列（要求与训练时的窗口长度一致），输出下一个事件的震级预测值。在实际运行时，需要维护一个循环缓冲队列，每来一条新地震目录记录就推进队列，丢弃最旧的事件，保持窗口长度恒定。这样可以做到每次新事件发生后自动更新下一个震级预测。

验证部署正确的方法很简单：取测试集的某一段连续序列，用 predictNextMagnitude 函数预测，和之前训练时的预测结果做对比。如果两次结果完全一致，说明模型加载、归一化参数、特征构造三个环节都正确；任何微小的偏差都说明某个环节与训练流程不一致。

部署时的另一个要点是异常处理。当窗口内包含缺失值（比如某次地震没有深度数据）时，可以直接用 fillmissing 补上上一事件的深度值。地震目录中这种缺失时有发生，不能因为一条数据就让整个预测流程崩溃。

本文还有配套的精品资源，点击获取