
简介该资源是一份面向环境科学专业学生、水务工程技术人员及研究人员的PDF文档聚焦区域级多井地下水位时空预测这一实际工程难题。其核心是融合图卷积网络与长短期记忆网络的GCN-LSTM模型通过构建观测井空间图结构结合空间自相似与属性自相似矩阵同步预测多口井的水位变化并引入温度、降雨等气象因素提升精度。资源包内共1个PDF文件大小约3.47MB完整呈现了模型原理、网络结构、数据集构建与实验分析包含成都城区56口井五年历史数据的验证过程。文档详细推导了GCN正向传播公式、空间与属性相似性计算以及编码器-解码器架构的设计思路并对比了仅考虑水位特征与融合多时间特征两种预测情景。目前已有436人学习适合希望掌握时空图神经网络在水文建模中应用、需要复现或借鉴多井同步预测方案的读者参考。1. 区域级地下水位预测为什么不能只靠单井 LSTM区域级多井地下水位预测本质上是一个带空间依赖的时空序列问题。单井 LSTM 能学好一条观测井的时间规律却学不到井与井之间的水力联系——上游抽水、下游水位跟着降这种空间传导单井模型完全看不见。GCN-LSTM 的思路是先用图卷积网络GCN在井网拓扑上做空间聚合再把聚合后的特征喂给 LSTM 捕捉时间依赖最后输出每口井未来若干天的水位。它解决的是「多井联合预测」而不是「逐井独立预测」适合水文监测站网、矿区沉降观测、灌区地下水位管理的从业者。下面这套流程我在实际项目里跑通过从建图、训练到落盘推理都有可复现的代码参数和踩坑点一并写清楚。2. GCN-LSTM 的建图逻辑与数据准备井网怎么变成邻接矩阵2.1 为什么用图结构表达井网空间关系地下水位在空间上不是孤立的。两口井距离越近、含水层连通性越好一口井的水位波动越容易传导到另一口井。GCN 的核心操作是邻接矩阵 A 与节点特征 X 的聚合每个节点井在每一层把邻居节点的特征加权求和权重由 A 决定。所以建图的质量直接决定空间建模的上限。常见的建图方式有三种基于地理距离的高斯核、基于水位序列相关性的动态图、以及两者融合的混合图。我一般用距离阈值加高斯核因为物理意义清晰、参数少、可解释。具体做法是计算所有井对的欧氏距离超过阈值 R 的置零阈值内的用 exp(-d²/σ²) 加权。R 和 σ 是两个必须调的参数R 取研究区井距中位数的 1.5 到 2 倍比较稳σ 取 R 的一半。注意邻接矩阵要做对称归一化否则度数大的节点特征会爆炸训练直接发散。2.2 数据格式与缺失值处理输入数据是一张长表每行是「井号 日期 水位」。宽表化之后得到形状为 (T, N) 的矩阵T 是时间步N 是井数。缺失值用线性插值加前后向填充不要用均值填充——水位是连续过程均值填充会破坏时间自相关LSTM 学出来的是错的。import numpy as np import pandas as pd # df: columns [well_id, date, water_level] df[date] pd.to_datetime(df[date]) pivot df.pivot_table(indexdate, columnswell_id, valueswater_level) pivot pivot.asfreq(D) # 统一为日尺度 pivot pivot.interpolate(methodlinear, limit_directionboth) pivot pivot.ffill().bfill() # 边界补齐 data pivot.values.astype(np.float32) # shape: (T, N) well_ids pivot.columns.tolist()这段代码做了三件事把长表转成 (T, N) 矩阵、按日重采样保证时间等间隔、用线性插值处理缺失。limit_directionboth保证首尾缺失也能补上。如果你的数据是月尺度把asfreq(D)改成asfreq(MS)即可但 LSTM 的窗口长度要相应调整。2.3 邻接矩阵构建代码from scipy.spatial.distance import cdist # coords: shape (N, 2), 每口井的经纬度或投影坐标 coords np.array([well_coords[w] for w in well_ids]) dist cdist(coords, coords, metriceuclidean) R np.median(dist[dist 0]) * 1.8 # 阈值井距中位数的1.8倍 sigma R / 2.0 A np.exp(-(dist ** 2) / (sigma ** 2)) A[dist R] 0 # 超阈值置零 np.fill_diagonal(A, 1.0) # 自环 # 对称归一化: D^{-1/2} A D^{-1/2} D np.diag(A.sum(axis1)) D_inv_sqrt np.linalg.inv(np.sqrt(D)) A_norm D_inv_sqrt A D_inv_sqrt A_norm A_norm.astype(np.float32)R控制图的稀疏度太小图会碎成多个不连通子图GCN 退化成单井模型太大所有井全连接空间信息被平均掉等于没建图。sigma控制权重衰减速度一般取 R 的一半。归一化那三行是标准做法别省。2.4 滑动窗口切样本def make_samples(data, A, input_len30, pred_len7): X, Y [], [] T data.shape[0] for t in range(T - input_len - pred_len 1): X.append(data[t:tinput_len]) # (input_len, N) Y.append(data[tinput_len:tinput_lenpred_len]) # (pred_len, N) return np.stack(X), np.stack(Y) X, Y make_samples(data, A_norm, input_len30, pred_len7) # X: (S, 30, N), Y: (S, 7, N)input_len30是回看 30 天pred_len7是预测未来 7 天。这两个值不是拍脑袋定的回看窗口至少要覆盖一个完整的水位响应周期7 天预测对应周级调度需求。如果你的区域有强季节性input_len 要拉到 90 以上。3. 模型搭建与训练GCN 和 LSTM 怎么串起来3.1 网络结构设计整体结构是输入 (batch, input_len, N) → 每个时间步做 GCN 空间聚合 → 得到 (batch, input_len, N, hidden) → 沿时间维送入 LSTM → 取最后隐状态 → 全连接输出 (batch, pred_len, N)。关键设计点GCN 是逐时间步共享权重的也就是说所有时间步用同一个邻接矩阵和同一套 GCN 参数。这样参数量小、不容易过拟合也符合「空间关系不随时间突变」的物理假设。import torch import torch.nn as nn 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, A): # x: (batch, N, in_dim), A: (N, N) x self.linear(x) x torch.einsum(nn, bnd - bnd, A, x) # 邻居聚合 return torch.relu(x) class GCNLSTM(nn.Module): def __init__(self, num_nodes, gcn_hidden32, lstm_hidden64, pred_len7): super().__init__() self.gcn1 GCNLayer(1, gcn_hidden) self.gcn2 GCNLayer(gcn_hidden, gcn_hidden) self.lstm nn.LSTM(gcn_hidden, lstm_hidden, batch_firstTrue) self.fc nn.Linear(lstm_hidden, pred_len) self.num_nodes num_nodes self.pred_len pred_len def forward(self, x, A): # x: (batch, input_len, N) b, t, n x.shape x x.permute(0, 2, 1).reshape(b * n, t, 1) # (b*n, t, 1) x x.reshape(b, n, t).permute(0, 2, 1) # (b, t, n) # 逐时间步 GCN outs [] for i in range(t): xi x[:, i, :].unsqueeze(-1) # (b, n, 1) xi self.gcn1(xi, A) xi self.gcn2(xi, A) # (b, n, gcn_hidden) outs.append(xi) h torch.stack(outs, dim1) # (b, t, n, gcn_hidden) h h.permute(0, 2, 1, 3).reshape(b * n, t, -1) lstm_out, _ self.lstm(h) last lstm_out[:, -1, :] # (b*n, lstm_hidden) out self.fc(last) # (b*n, pred_len) out out.reshape(b, n, self.pred_len).permute(0, 2, 1) return out # (b, pred_len, n)gcn_hidden32是空间特征维度lstm_hidden64是时间隐状态维度。这两个值在 N 小于 50 的井网上够用井数上百时 gcn_hidden 可以加到 64但要注意过拟合。einsum(nn, bnd - bnd, A, x)就是邻接矩阵乘特征等价于对每个节点做邻居加权求和。3.2 训练循环与损失函数device torch.device(cuda if torch.cuda.is_available() else cpu) model GCNLSTM(num_nodeslen(well_ids)).to(device) A_tensor torch.tensor(A_norm).to(device) optimizer torch.optim.Adam(model.parameters(), lr1e-3, weight_decay1e-5) criterion nn.MSELoss() X_t torch.tensor(X).to(device) # (S, 30, N) Y_t torch.tensor(Y).to(device) # (S, 7, N) dataset torch.utils.data.TensorDataset(X_t, Y_t) loader torch.utils.data.DataLoader(dataset, batch_size32, shuffleTrue) for epoch in range(100): model.train() total_loss 0 for xb, yb in loader: optimizer.zero_grad() pred model(xb, A_tensor) loss criterion(pred, yb) loss.backward() torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm5.0) optimizer.step() total_loss loss.item() if (epoch 1) % 10 0: print(fEpoch {epoch1}, Loss: {total_loss/len(loader):.6f})clip_grad_norm_那行是后悔药——LSTM 遇到水位突变时梯度容易炸不裁剪的话 loss 会突然变 NaN。weight_decay1e-5是轻量正则井数少的时候可以调到 1e-4。batch_size 取 32 是折中显存够可以加到 64。3.3 数据标准化与反标准化水位数值量纲差异大不标准化 LSTM 收敛很慢。用训练集的均值和标准差做 z-score验证集和测试集用同一套参数。train_end int(len(X) * 0.7) val_end int(len(X) * 0.85) mean X[:train_end].mean() std X[:train_end].std() 1e-8 X_norm (X - mean) / std Y_norm (Y - mean) / std # 预测后反标准化 pred_real pred * std mean提示mean 和 std 必须只用训练集算用了全量数据就是信息泄漏验证指标会虚高。4. 避坑与排查GCN-LSTM 训练中翻车的五个场景4.1 Loss 不下降反而震荡现象训练 loss 在前几个 epoch 下降之后开始上下大幅震荡验证 loss 持续升高。原因学习率太大或者邻接矩阵归一化没做对导致特征尺度不一致。GCN 聚合后节点特征方差会随度数变化没归一化时高度数节点输出值远大于低度数节点。解决先把 lr 降到 1e-4 试一轮确认 A_norm 每行和接近 1在 GCN 层后加 LayerNorm。4.2 预测值全部趋近于均值现象模型输出的未来 7 天水位几乎是常数和输入序列的波动完全无关。原因LSTM 隐状态维度太小或者 input_len 太短模型学不到有效时间模式退化成预测均值。另一个常见原因是损失函数被大量平稳井主导波动大的井贡献被淹没。解决把 lstm_hidden 从 64 加到 128input_len 从 30 加到 60对每口井的 loss 做加权权重取该井水位标准差的倒数。4.3 某些井预测误差特别大现象整体 MAE 看起来还行但个别井的预测误差是其他井的 5 到 10 倍。原因这些井在邻接矩阵里是孤立节点或弱连接节点GCN 拿不到有效空间信息等于只靠 LSTM 单井预测。也可能是这些井本身缺失值太多插值填充引入了虚假模式。解决检查 A_norm 中这些井的度数如果小于 2 就放宽 R对缺失率超过 30% 的井考虑在训练时降权或直接剔除。4.4 验证集 loss 远高于训练集现象训练 loss 降到 0.001验证 loss 停在 0.01 下不去。原因过拟合。井数少、样本少的时候 GCN-LSTM 参数量相对过剩。另外滑动窗口切样本时相邻样本高度重叠训练集和验证集如果随机划分会泄漏。解决按时间顺序划分不要随机打乱加 dropoutLSTM 层设 0.2weight_decay 加到 1e-4减少 gcn_hidden。4.5 推理时显存溢出现象训练时正常推理时 batch 一大就 OOM。原因GCNLSTM 的 forward 里对每个时间步循环做 GCN时间步长时中间激活值累积。input_len90 时显存占用是 input_len30 的三倍。解决推理时用torch.no_grad()把 batch_size 降到 16或者把逐时间步 GCN 改成先 reshape 再一次性矩阵乘减少中间变量。5. 进阶技巧用残差连接和动态图提升区域预测精度5.1 残差 GCN 缓解过平滑GCN 堆两层以上会出现过平滑——所有节点特征趋同空间区分度消失。加残差连接是最省事的解法每层 GCN 输出加上输入。class ResidualGCN(nn.Module): def __init__(self, dim): super().__init__() self.gcn GCNLayer(dim, dim) self.norm nn.LayerNorm(dim) def forward(self, x, A): return self.norm(x self.gcn(x, A))把原来 GCNLSTM 里的 gcn2 换成 ResidualGCN训练稳定性和预测精度都会有可见提升。LayerNorm 放在残差之后保证输出尺度一致。5.2 动态图让邻接矩阵随时间变化固定邻接矩阵假设空间关系不随时间变但实际中季节性抽水、灌溉周期会让井间相关性发生漂移。动态图的做法是用滑动窗口算每段时间的井间相关系数和距离图加权融合。def dynamic_adj(data_window, A_dist, alpha0.5): # data_window: (window_len, N) corr np.corrcoef(data_window.T) corr np.nan_to_num(corr, nan0.0) corr (corr 1) / 2 # 映射到 [0,1] np.fill_diagonal(corr, 1.0) A_dyn alpha * A_dist (1 - alpha) * corr D np.diag(A_dyn.sum(axis1)) D_inv_sqrt np.linalg.inv(np.sqrt(D 1e-8)) return (D_inv_sqrt A_dyn D_inv_sqrt).astype(np.float32)alpha控制距离图和相关图的权重0.5 是起点。窗口长度取 60 到 90 天太短相关系数噪声大太长反映不出动态变化。每个预测步用对应窗口的动态图推理时也要同步更新。5.3 评估指标与验证方法不要只看 MAE。区域级预测要同时看三个指标指标含义合格线参考MAE平均绝对误差小于水位日变幅的 20%RMSE均方根误差小于水位日变幅的 30%NSE纳什效率系数大于 0.75NSE 是水文领域最认的指标大于 0.75 算可用大于 0.85 算好。验证时按时间顺序留出最后 15% 做测试集不要随机抽。另外建议做一次「留一井交叉验证」每次拿掉一口井不参与训练看模型能不能靠邻居井预测它这能直接检验空间建模是否真的有效。我自己的习惯是每次调完参数先跑一遍留一井验证如果拿掉某口井后 NSE 掉到 0.5 以下说明这口井在图上太孤立得回头检查建图参数。这套流程跑顺之后区域级多井预测的精度比逐井 LSTM 通常能提升 15% 到 30%井网越密提升越明显。希望帮到你。本文还有配套的精品资源点击获取