ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

溶解氧预测进阶:从LSTM到EMD-LSTM信号分解实战指南

溶解氧预测进阶:从LSTM到EMD-LSTM信号分解实战指南 简介面向环境监测与水处理研究场景基于深度学习的溶解氧时间序列预测模型被整理为完整工程包涵盖数据预处理、信号分解、模型构建与异常检测等关键环节。包内共16个文件包括11个Python脚本和5个CSV数据文件脚本实现了EMD-LSTM、EEMD-LSTM等预测模型以及DBSCAN、孤立森林、LOF、OneClassSVM等异常检测算法CSV则提供了多个水域的实测溶解氧记录可供模型直接训练与验证。整体压缩包约925KB已有156人学习下载。对于具备一定Python基础、希望快速进入水质预测研究的环境科学学习者通过阅读源码和运行自带数据可复现完整建模流程深入理解LSTM系列网络在长时序依赖建模中的设计思路学习缺失值填充、数据归一化及异常样本剔除的标准处理方法为后续模型改进或实际工程应用提供可扩展的代码基底。1. 溶解氧预测为什么需要从 LSTM 走向信号分解水质监测里溶解氧DO浓度是最难做时序预测的指标之一。它的变化同时受水温、光照、藻类活动、复氧速率等多重因素影响曲线里既有昼夜周期性又有降雨后的突变段单纯用原始序列喂给 LSTM模型往往学到的是均值回归而不是真正的动态模式。项目里除了LSTM.py之外还放进了EMD-LSTM.py、EEMD-LSTM.py、EEMD_BP.py和eemd.py这个结构本身就是一条清晰的技术路线先用经验模态分解把非平稳的 DO 序列拆成本征模态函数IMF再对每个分量做深度学习建模最后叠加输出。这套组合拳对溶解氧这类强非线性、非平稳的水质参数往往比端到端的 LSTM 提升 10% 到 20% 的 RMSE 表现。适合正在做水质预测课题的学生以及需要在污水厂、地表水监测场景里落地预测算法的工程师。下文按项目文件顺序把这套流程从数据处理到组合模型复现一遍。2. 预处理与异常值截断分解之前先做数据体检2.1 缺失值填充与时间轴对齐打开PreProcessing.py能看到作者处理水质记录的基本套路。原始 CSV 文件里常见的问题有三种时间戳不连续、溶解氧字段含空值、同一采样点记录了多行重复数据。第一步通常是重采样和对齐把时间轴统一到固定频率比如把siwan_wuwan.csv和Water Quality Record.csv里的记录按小时重采样。常见做法是先用 pandas 把时间列解析成 datetime 索引再resample到目标频率缺失的观测值用插值补上。import pandas as pd import numpy as np df pd.read_csv(Water Quality Record.csv, parse_dates[time], index_coltime) df df.resample(H).mean() # 统一到小时频率 # 缺失值插值DO浓度通常平滑变化线性插值比前向填充更合理 df[do] df[do].interpolate(methodlinear, limit_directionboth) df[do] df[do].fillna(methodffill).fillna(methodbfill) # 删除 still 缺失的孤立点 df df.dropna(subset[do])逻辑说明resample(H).mean()把同一小时内的多点观测压缩成一个平均值interpolate用线性插值补齐中间缺口ffill/bfill处理序列首尾的空洞。这里有一个容易踩的坑如果原始数据里某个时间段完全无记录resample会生成长段 NaN此时线性插值会“无中生有”造出平滑曲线掩盖真实波动。我一般会加一个阈值判断连续缺失超过 3 个小时就直接删除该段而不是强行插值。2.2 四种异常检测算法的选型依据项目OutlierDetection目录里放了四个独立脚本DBSCAN.py、OneClassSVM.py、LOF.py、IF.py。这不是为了炫技而是因为溶解氧的“异常”定义并不唯一。传感器漂移会产生孤立尖峰暴雨径流会被检测为局部离群而水体翻塘事件则是持续整段的异常。四种方法在 sklearn 里都是一两行能调用的但适用场景差别很大算法核心思想适合的异常类型关键参数对 DO 场景的适用性LOF局部密度对比局部尖峰、孤立突变n_neighbors20适合识别传感器瞬时毛刺IF随机切割特征空间高维稀疏离群n_estimators100适合多变量联合异常DO温度pHOneClassSVM学习决策边界分布漂移nu0.05,gammascale适合整体分布偏移检测DBSCAN密度连通任意形状的离群簇eps0.5,min_samples5适合识别持续异常段from sklearn.neighbors import LocalOutlierFactor # DO 序列是一维的但可以和温度组成二维特征空间 X df[[do, temp]].values # contamination 是预估的异常比例水质数据一般取 1% - 5% lof LocalOutlierFactor(n_neighbors20, contamination0.03) outlier_mask lof.fit_predict(X) -1 # 不直接删除而是先标记人工抽检后再决定 df[outlier_flag] outlier_mask df.loc[outlier_mask, do] np.nan df[do] df[do].interpolate(methodcubic)参数说明n_neighbors决定局部邻域的大小水质数据采样频率高时取 20 到 50 比较稳定contamination0.03意思是先验认为 3% 的数据是异常的这个值需要根据实际传感器运维记录调整调太小会漏掉尖峰调太大把真实低氧事件当异常删掉。对 DO 序列的多维检测推荐把温度、pH、电导率一起作为特征异常检测的目标是排除“传感器读数异常”而不是排除“水质本身异常”的真实低氧事件——两者存在本质区别。2.3 归一化与训练集划分异常值处理后PreProcessing.py里通常还会做滚动窗口划分。时间序列不能随机打乱必须按时间顺序切分。常见做法是前 70% 做训练中间 15% 做验证最后 15% 做测试这个比例对水质日尺度数据比较合理。归一化方法我建议对 DO 序列用MinMaxScaler对温度等辅助特征用StandardScaler原因在于后续 EMD 分解要求输入值域不过度扭曲MinMax 能把 DO 映射到 [0,1] 区间减少 LSTM 的激活函数饱和问题而温度特征分布本身接近正态标准化更合适。from sklearn.preprocessing import MinMaxScaler, StandardScaler from sklearn.model_selection import train_test_split scaler_do MinMaxScaler(feature_range(0, 1)) scaler_temp StandardScaler() do_scaled scaler_do.fit_transform(df[[do]]) temp_scaled scaler_temp.fit_transform(df[[temp]]) features np.hstack([do_scaled, temp_scaled]) # 时间序列切分不能用 train_test_split 随机打乱 train_end int(len(features) * 0.7) val_end int(len(features) * 0.85) train, val, test features[:train_end], features[train_end:val_end], features[val_end:]逻辑说明train_test_split默认是随机抽样对时间序列是灾难性的模型会在训练时看到未来数据验证集失去意义。手写切分保证时序完整性。MinMaxScaler的feature_range(0,1)和 LSTM 的 tanh 输出域匹配。3. LSTM 基线与训练配置先让单模型跑通3.1 滑动窗口构造监督样本LSTM.py实现的是标准的 LSTM 序列预测。在 DO 预测场景里输入特征是过去 24 小时到 72 小时的观测序列输出是未来 1 小时或 6 小时的 DO 浓度。滑窗构造是这类任务的核心环节直接把一维序列转换成[样本数, 时间步长, 特征数]的三维张量。def create_sequences(data, seq_len24, pred_len1): X, y [], [] for i in range(len(data) - seq_len - pred_len 1): X.append(data[i:iseq_len, :]) y.append(data[iseq_len:iseq_lenpred_len, 0]) # 只预测 DO return np.array(X), np.array(y) seq_len 24 # 用过去24小时的数据 pred_len 6 # 预测未来6小时 X_train, y_train create_sequences(train, seq_len, pred_len) X_val, y_val create_sequences(val, seq_len, pred_len)逻辑说明pred_len6意味着要预测未来 6 个时点这比单步预测更有实际意义——水务调度需要提前几小时知道溶解氧会不会跌破阈值。但多步预测有两种实现方式一种是这里展示的直接多输出另一种是递归预测把上一步输出拼回输入。直接多输出训练更稳定递归预测误差会随时间累积。对 DO 浓度这种中等平稳性的序列直接多输出更合适。3.2 网络结构与训练超参LSTM 层数不宜过深水质时间序列样本量通常只有几千到几万条叠三层以上极容易过拟合。项目里Sever_train_version.py应该是整合了调参后的完整版常见配置是单层或双层 LSTM隐藏单元 32 到 64dropout 0.2配合早停机制。训练参数推荐如下from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense, Dropout from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau model Sequential([ LSTM(64, return_sequencesTrue, input_shape(seq_len, X_train.shape[2])), Dropout(0.2), LSTM(32, return_sequencesFalse), Dropout(0.2), Dense(pred_len) ]) model.compile(optimizeradam, lossmse, metrics[mae]) callbacks [ EarlyStopping(monitorval_loss, patience10, restore_best_weightsTrue), ReduceLROnPlateau(monitorval_loss, factor0.5, patience5, min_lr1e-5) ] history model.fit( X_train, y_train, validation_data(X_val, y_val), epochs100, batch_size64, callbackscallbacks, verbose0 )参数说明EarlyStopping的patience10表示验证损失连续 10 个 epoch 不下降就停restore_best_weightsTrue回滚到最优权重。ReduceLROnPlateau在验证集停滞时把学习率减半避免震荡。batch_size 取 64 是内存和梯度稳定性的折中样本量小时 batch_size 设 16 到 32 能提升泛化。lossmse的平方惩罚会让模型更关注偏离大的点这对 DO 预测是双刃剑——可能过度拟合尖峰噪声实际项目中我倾向用huber_loss降低离群点权重。3.3 损失函数与评估指标的坑项目摘要里提到了 MSE 和 RMSE但实际评估 DO 预测模型时RMSE 单指标会掩盖系统性偏差。比如模型预测整体偏低但趋势正确RMSE 可能很小但对“溶解氧低于 2 mg/L 的缺氧事件”完全无法预警。我一般会额外计算 MAE 和 R²同时把测试集按 DO 阈值分段统计误差from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score y_pred model.predict(X_test) mae mean_absolute_error(y_test, y_pred) rmse np.sqrt(mean_squared_error(y_test, y_pred)) r2 r2_score(y_test, y_pred) # 按溶解氧等级拆分误差 low_oxygen y_test 2 print(f缺氧段 MAE: {mean_absolute_error(y_test[low_oxygen], y_pred[low_oxygen]):.3f}) print(f正常段 MAE: {mean_absolute_error(y_test[~low_oxygen], y_pred[~low_oxygen]):.3f})这里的y_test 2是缺氧阈值生活饮用水水源地标准一般要求 DO 不低于 3 mg/L养殖水体则各有差异。分段评估能直接看出模型在关键区间的表现这是论文里不会写但工程上很重要的一步。4. EMD/EEMD 分解与组合模型的可复现实现4.1 EMD 分解的基本逻辑与模态混叠eemd.py是整套代码的核心也是和前两章内容连接的关键。经验模态分解EMD把原始 DO 序列分解为若干个 IMF 和一个残差项每个 IMF 代表不同频率的振动模态。高频 IMF 对应短时波动和噪声低频 IMF 和残差对应昼夜循环和季节趋势。这样 LSTM 只需要对不同频率的平稳分量建模难度大幅降低。但标准 EMD 有一个严重缺陷就是模态混叠——相近频率的分量被分散到不同 IMF 里导致分解结果不稳定。集合经验模态分解EEMD通过在原序列中加入白噪声再进行多次平均来解决这个问题。项目里同时出现emd.py和eemd.py说明作者做了对比实验。# 基于 PyEMD 库实现 EEMD这是最常见的解法 from PyEMD import EEMD eemd EEMD() eemd.noise_seed(42) # 这两个参数是 EEMD 的核心控制项 eemd_noise 0.2 # 添加白噪声的幅度占原始信号标准差比例 eemd_trials 100 # 集成次数 imfs eemd.eemd(do_series, noise_widtheemd_noise, ensemble_sizeeemd_trials)逻辑说明noise_width0.2表示加入的白噪声标准差为原始信号的 20%太小无法有效改变极值分布太大会把噪声当成真实信号ensemble_size100是集成平均次数增大能有效抑制噪声残留但计算开销线性增长。对 DO 这类日周期明显的数据我一般先看分解出的 IMF 数量是否在 6 到 9 个之间如果低于 5 个说明序列本身不够复杂没必要用 EEMD直接 LSTM 就能建模。4.2 EEMD-LSTM 的多分量建模与重构EEMD_LSTM.py的实现策略通常是对每个 IMF 单独训练一个 LSTM然后叠加预测结果。这是经典做法但工程上有个效率问题假设分解出 8 个 IMF就要训练 8 个 LSTM训练时间和调参工作量线性放大。我一般会在代码里加一个判断先对各个 IMF 做样本熵或方差贡献率计算只对贡献率超过 1% 的分量建模其余合并到残差项。另一种做法是把所有 IMF 拼成多维输入一个 LSTM 多输出但实测效果通常不如独立建模好——不同 IMF 的时间尺度差异太大共享网络难以同时拟合高频和低频特征。from PyEMD import EEMD from tensorflow.keras.models import Sequential from tensorflow.keras.layers import LSTM, Dense def eemd_lstm_predict(series, n_imfsNone): eemd EEMD() imfs eemd.eemd(series, noise_width0.2, ensemble_size100) predictions np.zeros(len(series)) imf_preds [] for i, imf in enumerate(imfs): # 构造滑窗样本 X, y create_sequences(imf.reshape(-1, 1), seq_len24, pred_len1) split int(len(X) * 0.8) model Sequential([ LSTM(32, input_shape(X.shape[1], 1)), Dense(1) ]) model.compile(optimizeradam, lossmse) model.fit(X[:split], y[:split], epochs30, batch_size32, validation_data(X[split:], y[split:]), verbose0) # 对末端做多步预测 last_seq imf[-24:].reshape(1, 24, 1) pred model.predict(last_seq, verbose0)[0, 0] imf_preds.append(pred) return np.sum(imf_preds), imf_preds, imfs逻辑说明这里对每个 IMF 用相同的seq_len24是简化处理实际调参时高频 IMF 用短窗口更灵敏低频分量用长窗口更能捕捉周期。imf_preds列表保存每个分量的单步预测值最后np.sum还原为完整预测。注意create_sequences需要复用第 3 章的函数保证滑窗逻辑一致。4.3 BP 与 LSTM 的对比价值EEMD_BP.py存在的意义不只是换个模型跑个对比图它是判断“分解贡献了多少、深度学习贡献了多少”的重要对照。BP 网络结构简单、没有时序记忆能力如果 EEMD-BP 能达到和 EEMD-LSTM 相近的精度说明分解后的 IMF 已经足够平稳LSTM 的长期记忆能力没有发挥出来此时可以简化模型结构甚至改用更轻量的回归算法反之如果 EEMD-LSTM 显著优于 EEMD-BP说明 IMF 内部仍存在时序依赖LSTM 的循环结构真正起了作用。from sklearn.neural_network import MLPRegressor from sklearn.preprocessing import StandardScaler def eemd_bp_predict(series): eemd EEMD() imfs eemd.eemd(series, noise_width0.2, ensemble_size100) imf_preds [] for imf in imfs: # BP 用滑窗拉平的特征没有时序结构 X, y create_sequences(imf.reshape(-1, 1), seq_len24, pred_len1) X_flat X.reshape(X.shape[0], -1) split int(len(X) * 0.8) scaler StandardScaler() X_scaled scaler.fit_transform(X_flat[:split]) X_val_scaled scaler.transform(X_flat[split:]) mlp MLPRegressor(hidden_layer_sizes(64, 32), max_iter200, learning_rate_init0.001, early_stoppingTrue) mlp.fit(X_scaled, y[:split]) # 用训练集均值填充验证集 pred mlp.predict(X_val_scaled[-1:])[0] imf_preds.append(pred) return np.sum(imf_preds)参数说明hidden_layer_sizes(64, 32)是两层隐藏层learning_rate_init0.001是 Adam 优化器的默认起点。MLP 的early_stoppingTrue同样监控验证损失防止过拟合。这里把seq_len24的窗口拉平成 24 维特征等于放弃了时间顺序信息这正好是 BP 和 LSTM 的核心差异。5. 用残差分布回溯验证分解带来的真实增益第 4 章的对比实验完成后还有一个容易被忽略的验证步骤检查每个 IMF 的预测误差是否服从零均值白噪声。如果某个 IMF 的残差序列仍然存在明显的趋势或自相关说明该分量的分解不充分直接加总会导致系统偏差累积。from statsmodels.stats.diagnostic import acorr_ljungbox def check_imf_residual(residuals, name): # Ljung-Box 检验残差是否为白噪声 lb_test acorr_ljungbox(residuals, lags10, return_dfTrue) p_value lb_test[lb_pvalue].iloc[-1] if p_value 0.05: print(f{name}: 残差为白噪声分解充分 (p{p_value:.3f})) else: print(f{name}: 残差仍有自相关需调整分解参数 (p{p_value:.3f}))# 对第 4 章的 imf_preds 结果做残差检查 # residuals 每个 IMF 实际值 - 预测值 for i, res in enumerate(residuals_list): check_imf_residual(res, fIMF{i1})逻辑说明acorr_ljungbox检验残差序列在滞后 10 阶内是否存在自相关p 值大于 0.05 说明残差是白噪声预测已经捕捉到了该分量的全部可学习模式。这个检验能定位“分解不够深”还是“LSTM 没学好”——如果高频 IMF 残差白噪声但低频 IMF 有自相关说明低频分量的滑窗长度太短需要把seq_len从 24 增加到 48 或 72。EEMD 的noise_width和ensemble_size两个参数也值得在验证后回看。我发现对 DO 序列noise_width0.2在大多数情况下是稳健的但如果分解出的第一个 IMF 方差占比异常高说明噪声幅值偏大如果频率混叠现象明显则先提高ensemble_size到 200代价是运行时间翻倍。这个回溯流程跑完模型才算真正达到可交付状态而不仅仅是测试集上 RMSE 好看。实际部署时把训练好的分解参数和网络权重固化到配置文件里新来的数据直接走同一套预处理、分解、预测管道整个溶解氧时间序列预测系统就可以稳定运行了。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进