多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

区域地下水位预测实战:GCN-LSTM 模型构建与避坑指南

区域地下水位预测实战:GCN-LSTM 模型构建与避坑指南 简介这份资源面向环境科学专业学生、水务工程技术人员及相关领域研究人员聚焦区域级多井地下水位时空预测这一实际工程难题。其核心是融合图卷积网络与长短期记忆网络的GCN-LSTM模型先以成都城区56口观测井构建空间图结构用邻接矩阵刻画井间空间距离与属性相似性再由GCN提取空间关联、LSTM挖掘时间依赖实现多井水位同步预测并引入温度、降雨等气象因素提升精度。资源包共1个PDF文件约3.47MB完整呈现模型原理、网络结构、数据集构建与实验验证过程包含GCN正向传播公式、空间与属性自相似矩阵推导、编码-解码器结构及五年1826天观测数据的训练测试划分方案。已有439人学习适合希望掌握时空图神经网络在水文建模中落地思路、复现多井预测实验的读者参考。1. 区域级地下水位预测为什么不能只靠单井时序模型地下水位预测这件事单井做时序建模看起来很简单一口井的历史水位丢进 LSTM滚动预测未来 7 天误差往往还不错。但一旦把范围放大到区域级——比如一个灌区、一个盆地、几十口甚至上百口监测井——单井模型就会集体翻车。原因不复杂地下水位不是孤立的相邻井之间存在明显的水力联系上游抽水会传导到下游补给区的变化会延迟影响排泄区。你只喂一口井的历史数据模型学不到这种空间耦合关系预测出来的曲线在丰水期还能看一到枯水期或集中开采期就完全失真。GCN-LSTM 就是冲着这个痛点来的。GCN图卷积网络负责在井与井之间传播空间信息LSTM 负责在时间维度上捕捉水位变化的惯性两者串起来才能同时吃到“空间相关性”和“时间依赖性”。这套方案适合做区域水资源评价的从业者、做地下水数值模拟想加一个数据驱动对照组的工程师以及手头有几十口井的监测数据、想快速搭一个可复现预测流程的人。它不替代 MODFLOW 这类物理模型但作为快速预警和缺测插补的工具性价比很高。2. 把监测井变成图GCN-LSTM 的建模逻辑与数据准备2.1 为什么用图结构而不是网格区域地下水位数据天然适合图结构。每口井是一个节点井与井之间的水力联系是边。相比把区域切成规则网格再卷积图结构不需要插值不会在无井区域引入虚假信息而且能直接处理不规则分布的监测网。常见做法是用井间距离加地质连通性来构造邻接矩阵距离越近、含水层连通性越好边权越大。如果只有距离信息用高斯核或反距离权重也能跑但预测精度会打折扣尤其是在隔水断层两侧的井距离近但水力联系弱不加约束容易带偏。GCN 的核心操作是聚合邻居信息。每一层图卷积节点会把自己的特征和邻居的特征按边权加权求和再过一个线性变换。堆两层就够了层数多了会过平滑所有井的特征趋同反而丢失局部差异。LSTM 接在 GCN 后面把每个时间步的图卷积输出当作输入序列学习时间上的演变规律。训练时用历史 N 天预测未来 M 天N 一般取 30 到 90M 取 7 到 30具体看区域水位响应快慢。2.2 数据清洗与缺失值处理地下水位监测数据最大的问题是缺测和异常。常见做法是先按井号和时间排序把明显超出物理范围的值比如水位突然跳变几十米标记为异常再用线性插值或 KNN 插补。注意不要用全局均值填充那会抹掉井间差异。对于连续缺测超过 7 天的井段建议直接剔除该时段不要硬插。import pandas as pd import numpy as np # 假设原始数据列well_id, date, water_level df pd.read_csv(groundwater_raw.csv, parse_dates[date]) df df.sort_values([well_id, date]) # 按井分组标记异常值超过该井均值 ±3 倍标准差 def mark_outliers(group): mean, std group[water_level].mean(), group[water_level].std() group[is_outlier] np.abs(group[water_level] - mean) 3 * std return group df df.groupby(well_id, group_keysFalse).apply(mark_outliers) # 异常值置为 NaN再线性插值 df.loc[df[is_outlier], water_level] np.nan df[water_level] df.groupby(well_id)[water_level].transform( lambda s: s.interpolate(methodlinear, limit7) ) # 剔除仍为 NaN 的行连续缺测超过 7 天 df df.dropna(subset[water_level])这段代码的逻辑是先按井计算统计量避免不同井水位基准不同导致误判异常值置 NaN 后限制插值长度防止长段缺测被过度平滑。参数limit7可以根据监测频率调整日监测数据用 7 比较稳小时数据可以放宽到 24。2.3 构造邻接矩阵与归一化邻接矩阵的质量直接决定 GCN 能不能学到有用的空间关系。如果只有井位坐标用反距离权重构造from scipy.spatial.distance import cdist coords df.groupby(well_id)[[lon, lat]].first().values dist cdist(coords, coords, metriceuclidean) # 反距离权重避免除零 sigma np.median(dist[dist 0]) adj np.exp(-(dist ** 2) / (2 * sigma ** 2)) np.fill_diagonal(adj, 0) # 对称归一化D^{-1/2} A D^{-1/2} degree adj.sum(axis1) d_inv_sqrt np.diag(1.0 / np.sqrt(degree 1e-8)) adj_norm d_inv_sqrt adj d_inv_sqrtsigma取距离中位数是经验做法能让权重衰减速度匹配监测网密度。如果井很密sigma 小一点井稀疏sigma 大一点。归一化是为了防止度数大的节点在聚合时数值爆炸这是 GCN 的标准操作不能省。3. 用 PyTorch 搭一个可复现的 GCN-LSTM 训练流程3.1 模型定义图卷积层与 LSTM 的拼接方式GCN-LSTM 有两种常见拼接一种是先 GCN 再 LSTM每个时间步做一次图卷积把结果按时间排成序列喂给 LSTM另一种是 GCN 和 LSTM 并行最后融合。前者更适合地下水位这种空间关系相对稳定、时间演变主导的场景。下面用 PyTorch 实现第一种。import torch import torch.nn as nn import torch.nn.functional as F class GCNLayer(nn.Module): def __init__(self, in_dim, out_dim): super().__init__() self.linear nn.Linear(in_dim, out_dim) def forward(self, x, adj): # x: (batch, num_nodes, in_dim), adj: (num_nodes, num_nodes) support self.linear(x) out torch.einsum(nm,bmf-bnf, adj, support) return F.relu(out) class GCNLSTM(nn.Module): def __init__(self, num_nodes, in_dim, gcn_hidden, lstm_hidden, out_dim): super().__init__() self.gcn1 GCNLayer(in_dim, gcn_hidden) self.gcn2 GCNLayer(gcn_hidden, gcn_hidden) self.lstm nn.LSTM( input_sizegcn_hidden * num_nodes, hidden_sizelstm_hidden, num_layers2, batch_firstTrue, dropout0.2 ) self.fc nn.Linear(lstm_hidden, out_dim * num_nodes) self.num_nodes num_nodes self.out_dim out_dim def forward(self, x, adj): # x: (batch, seq_len, num_nodes, in_dim) batch, seq_len, num_nodes, in_dim x.shape gcn_outs [] for t in range(seq_len): h self.gcn1(x[:, t], adj) h self.gcn2(h, adj) gcn_outs.append(h) # (batch, seq_len, num_nodes * gcn_hidden) gcn_seq torch.stack(gcn_outs, dim1).reshape(batch, seq_len, -1) lstm_out, _ self.lstm(gcn_seq) pred self.fc(lstm_out[:, -1]) # 取最后时间步 return pred.reshape(batch, self.out_dim, self.num_nodes)关键点torch.einsum(nm,bmf-bnf, adj, support)实现的是邻接矩阵与节点特征的批量乘法比循环节点快很多。LSTM 输入维度是gcn_hidden * num_nodes因为每个时间步把所有节点的图卷积输出展平成一个向量。输出层再 reshape 回(batch, out_dim, num_nodes)表示每个节点预测未来 out_dim 个时间步。dropout0.2是防过拟合的常规设置数据量少可以调到 0.3。3.2 滑动窗口构造训练样本时间序列预测要把连续数据切成滑动窗口。假设用 60 天历史预测未来 7 天def create_sequences(data, input_len60, pred_len7): # data: (total_days, num_nodes, features) xs, ys [], [] for i in range(len(data) - input_len - pred_len 1): xs.append(data[i:i input_len]) ys.append(data[i input_len:i input_len pred_len, :, 0]) # 只预测水位 return np.array(xs), np.array(ys) # 假设 data 形状 (T, N, F)F 包含水位、降雨、温度等 X, Y create_sequences(data, input_len60, pred_len7) # 按时间划分训练/验证/测试不要随机打乱 split int(len(X) * 0.7) X_train, Y_train X[:split], Y[:split] X_val, Y_val X[split:int(len(X)*0.85)], Y[split:int(len(X)*0.85)] X_test, Y_test X[int(len(X)*0.85):], Y[int(len(X)*0.85):]注意划分必须按时间顺序随机打乱会导致未来信息泄漏验证集误差虚低上线就翻车。input_len和pred_len要根据区域水位响应时间调响应慢的区域 input_len 可以拉到 90响应快的 30 就够。3.3 训练循环与早停策略device torch.device(cuda if torch.cuda.is_available() else cpu) model GCNLSTM(num_nodesN, in_dimF, gcn_hidden32, lstm_hidden64, out_dim7).to(device) adj_tensor torch.tensor(adj_norm, dtypetorch.float32).to(device) optimizer torch.optim.Adam(model.parameters(), lr1e-3, weight_decay1e-4) criterion nn.MSELoss() best_val float(inf) patience, counter 15, 0 for epoch in range(200): model.train() train_loss 0 for i in range(0, len(X_train), 32): xb torch.tensor(X_train[i:i32], dtypetorch.float32).to(device) yb torch.tensor(Y_train[i:i32], dtypetorch.float32).to(device) optimizer.zero_grad() pred model(xb, adj_tensor) loss criterion(pred, yb) loss.backward() torch.nn.utils.clip_grad_norm_(model.parameters(), 5.0) optimizer.step() train_loss loss.item() model.eval() with torch.no_grad(): xv torch.tensor(X_val, dtypetorch.float32).to(device) yv torch.tensor(Y_val, dtypetorch.float32).to(device) val_loss criterion(model(xv, adj_tensor), yv).item() if val_loss best_val: best_val val_loss torch.save(model.state_dict(), best_gcn_lstm.pth) counter 0 else: counter 1 if counter patience: print(fEarly stop at epoch {epoch}) break梯度裁剪clip_grad_norm_(..., 5.0)是 LSTM 训练必备防止梯度爆炸。早停 patience 设 15 到 20太小容易停在局部最优太大浪费算力。学习率 1e-3 配 Adam 是稳妥起点如果 loss 震荡明显降到 5e-4。4. 避坑与排查GCN-LSTM 在地下水位场景的 5 个血泪教训4.1 邻接矩阵太稠密导致过平滑现象训练 loss 正常下降但验证集上所有井的预测值趋同空间差异消失。原因邻接矩阵没有做稀疏化每口井都跟所有其他井有不可忽略的边权GCN 聚合后节点特征被平均化。解决对邻接矩阵做阈值截断只保留每个节点的 top-k 邻居k 取 5 到 10。或者用距离阈值超过一定距离的边直接置零。代码里加一行adj[adj np.percentile(adj, 70)] 0再归一化。4.2 输入特征里混入未来信息现象验证集误差极低测试集误差突然翻倍。原因构造样本时用了全局归一化均值和方差包含了测试集信息或者滑动窗口切分时不小心让训练集包含了未来时间点。解决归一化统计量只用训练集计算然后应用到验证和测试集。滑动窗口严格按时间顺序切训练集的时间范围必须完全早于验证集。4.3 LSTM 层数过多导致训练不收敛现象loss 震荡验证集误差始终高于训练集很多。原因LSTM 堆了 3 层以上参数量过大而地下水位数据样本量通常只有几千条严重过拟合。解决LSTM 层数控制在 1 到 2 层hidden_size 从 64 起步不要一上来就 256。加 dropout 和 weight_decay如果还不行就减层。4.4 缺失值插补引入虚假趋势现象某口井连续缺测一个月插补后模型在这口井上预测异常平滑但实际恢复监测后误差很大。原因线性插值在长段缺测上会造出一条直线模型学到了这个虚假模式。解决连续缺测超过 7 天的段直接剔除不要插补。如果必须保留该井用该井历史同期的均值填充并加一个缺失标记特征让模型自己学。4.5 评估指标只看 RMSE 忽略峰值现象RMSE 看起来不错但枯水期或集中开采期的水位骤降完全没预测到。原因RMSE 对整体拟合敏感对峰值不敏感。地下水位预测最关键的恰恰是极值段。解决加 NSE纳什效率系数和峰值误差指标。NSE 大于 0.7 才算可用峰值误差单独看。训练时可以对极值段加权或者用 Huber loss 替代 MSE。5. 让预测结果真正可用滚动更新与不确定性输出模型训完只是第一步落地时最有用的是滚动预测和不确定性区间。地下水位预测如果只给一条曲线决策者没法判断风险。我一般会在 GCN-LSTM 外面套一层滚动更新每天用最新监测数据替换窗口里最旧的一天重新推理这样模型能跟上实时变化。实现上不需要重新训练只做前向传播成本很低。def rolling_predict(model, latest_window, adj, steps7): # latest_window: (1, input_len, N, F) model.eval() preds [] window latest_window.copy() with torch.no_grad(): for _ in range(steps): x torch.tensor(window, dtypetorch.float32).to(device) pred model(x, adj).cpu().numpy() # (1, 1, N) preds.append(pred[0, 0]) # 把预测值拼到窗口末尾去掉最早一天 new_day window[:, -1].copy() new_day[:, :, 0] pred[0, 0] # 水位列替换为预测 window np.concatenate([window[:, 1:], new_day[:, None]], axis1) return np.array(preds) # (steps, N)不确定性输出用 MC Dropout推理时保持 dropout 开启跑 50 次取均值和分位数。这样能给出 90% 置信区间比单点预测有用得多。注意 MC Dropout 的 dropout 率要和训练时一致否则区间偏窄。def mc_dropout_predict(model, x, adj, n_samples50): model.train() # 保持 dropout 开启 preds [] with torch.no_grad(): for _ in range(n_samples): preds.append(model(x, adj).cpu().numpy()) preds np.array(preds) mean preds.mean(axis0) lower np.percentile(preds, 5, axis0) upper np.percentile(preds, 95, axis0) return mean, lower, upper这套组合跑下来区域级预测的 NSE 通常能到 0.75 以上比单井 LSTM 高 0.1 到 0.15。但别指望它替代物理模型它的价值在于快和可复现。我自己的习惯是每次拿到新区域的数据先画井位图和邻接矩阵热力图确认空间结构合理再训模型。这一步花十分钟能省掉后面几天的排查。希望帮到你。本文还有配套的精品资源点击获取
返回列表