GCN-LSTM多井地下水位预测：从模型拆解到工程避坑实战

简介该资源是一份面向环境科学专业学生、水务工程技术人员及研究人员的PDF文档聚焦区域级多井地下水位时空预测难题。其核心是融合图卷积网络与长短期记忆网络的GCN-LSTM模型通过构建观测井空间图结构结合空间自相似与属性自相似矩阵同步预测多口井的水位变化并引入温度、降雨等气象因素提升精度。资源包内仅含1个PDF文件大小约3.47MB完整呈现了模型原理、网络结构、数据集构建及成都56口井五年实测数据的实验验证过程。已有439人学习下载。读者可从中获取从空间特征提取到时间序列建模的完整技术路线理解编码-解码器结构、邻接矩阵设计及多情景对比实验方法为城市水务管理与防灾减灾决策提供可复现的算法参考与理论支撑。1. 从单井到区域为什么GCN-LSTM成了多井地下水位预测的刚需如果你手头只有几口井的历史水位数据用ARIMA或者单变量LSTM跑一跑短期趋势基本能看。但一旦把范围拉到区域级——比如一个城市几十口甚至上百口观测井——单井模型就露怯了。原因不复杂地下水位不是孤立的相邻井之间存在明显的水力联系上游抽水会传导到下游降雨入渗在空间上也有扩散效应。传统时序模型把每口井当成独立序列处理等于主动扔掉了空间关联信息预测精度自然上不去。GCN-LSTM要解决的就是这个痛点。它先用图卷积网络GCN把观测井之间的空间拓扑关系编码进节点特征再用长短期记忆网络LSTM捕捉时间维度的演变规律最终一次性输出多口井的水位预测值。这套思路在区域级水务管理、地下水资源评估、防灾减灾辅助决策等场景里都有直接落地价值。适合环境科学、水文地质、水务工程方向的研究人员和工程技术人员尤其是手头有多年多井观测数据、想做时空联合建模的从业者。这份资源把模型结构、数据集构建、两种预测情景的对比实验都写清楚了复现路径比较完整。2. GCN-LSTM模型拆解从邻接矩阵到时空特征融合2.1 图结构怎么建空间自相似矩阵与属性自相似矩阵GCN的核心输入是图结构具体到地下水位预测场景节点就是观测井边就是井与井之间的关联强度。这个关联强度怎么算直接决定了GCN能提取到什么空间特征。原文给了两个矩阵空间自相似矩阵W和属性自相似矩阵C最终邻接矩阵A由两者共同决定。空间自相似矩阵W用的是欧式距离平方的反比。公式很直白W_ij 1/d²当i≠j时W_ij 0当ij时。这里有个工程上的简化假设——井与井之间的海拔差远小于水平距离所以连线与水平面的夹角忽略不计。这个假设在平原或盆地城市基本成立但如果你的研究区域地形起伏剧烈比如山区这个简化就会引入误差需要把高程差纳入距离计算。属性自相似矩阵C用的是余弦相似度。每口井有四个属性经度、纬度、孔口标高、井深。把这四个属性组成一个向量两口井之间的余弦值就是属性相似度。这个设计的逻辑是即使两口井空间距离不远但如果孔口标高和井深差异很大它们的水位关联性也可能不强。余弦相似度能捕捉这种属性层面的相似性。import numpy as np from sklearn.metrics.pairwise import cosine_similarity def build_spatial_matrix(coords): coords: shape (N, 2), 每口井的经度、纬度 返回空间自相似矩阵 W, shape (N, N) N coords.shape[0] W np.zeros((N, N)) for i in range(N): for j in range(N): if i ! j: # 欧式距离平方 d2 np.sum((coords[i] - coords[j]) ** 2) W[i, j] 1.0 / d2 if d2 0 else 0 return W def build_attribute_matrix(attrs): attrs: shape (N, 4), 每口井的经度、纬度、孔口标高、井深 返回属性自相似矩阵 C, shape (N, N) # 余弦相似度直接调用sklearn C cosine_similarity(attrs) # 对角线置零避免自环 np.fill_diagonal(C, 0) return C def build_adjacency(W, C, alpha0.5): 融合空间和属性相似度alpha控制两者权重 A alpha * W_norm (1 - alpha) * C_norm W_norm W / (W.sum(axis1, keepdimsTrue) 1e-8) C_norm C / (C.sum(axis1, keepdimsTrue) 1e-8) A alpha * W_norm (1 - alpha) * C_norm return A上面这段代码里build_spatial_matrix用的是距离平方反比注意当两口井坐标完全相同时会除零所以加了d2 0的判断。build_attribute_matrix直接调cosine_similarity但要把对角线置零否则每口井和自己相似度为1会形成自环影响GCN的传播效果。build_adjacency里的alpha是个可调参数默认0.5表示空间和属性权重各半。如果你的场景里空间距离影响更大可以把alpha调到0.7以上如果属性差异更关键就降到0.3左右。这个参数没有理论最优值需要根据验证集表现来调。2.2 GCN层的前向传播两层网络与激活函数选择原文给了一个两层的GCN网络第一层用ReLU第二层用Softmax。前向传播公式是Z softmax( · ReLU( · X · W⁰) · W¹)其中Â是归一化后的邻接矩阵X是节点特征矩阵W⁰和W¹是可学习的权重矩阵。这里有个细节值得注意第二层用Softmax通常是为了做节点分类但地下水位预测是回归任务输出应该是连续值。所以实际实现时第二层一般不加Softmax直接线性输出损失函数用MSE而不是交叉熵。原文的公式可能是从分类任务迁移过来的通用形式复现时要注意这个区别。Â的归一化方式原文写的是D⁻¹ÃD⁻¹这是对称归一化的一种变体。标准的GCN归一化是D⁻¹/²ÃD⁻¹/²两者都能用但对称归一化在训练时更稳定。如果你发现训练loss震荡厉害可以换成对称归一化试试。import torch import torch.nn as nn import torch.nn.functional as F class GCNLayer(nn.Module): def __init__(self, in_features, out_features, activationNone): super().__init__() self.linear nn.Linear(in_features, out_features) self.activation activation def forward(self, x, adj_norm): # x: (N, in_features), adj_norm: (N, N) support self.linear(x) # X * W output torch.mm(adj_norm, support) #  * (X * W) if self.activation is not None: output self.activation(output) return output class GCNEncoder(nn.Module): def __init__(self, in_dim, hidden_dim, out_dim): super().__init__() self.gcn1 GCNLayer(in_dim, hidden_dim, activationF.relu) self.gcn2 GCNLayer(hidden_dim, out_dim, activationNone) # 回归任务不加softmax def forward(self, x, adj_norm): h self.gcn1(x, adj_norm) h self.gcn2(h, adj_norm) return hGCNLayer里把线性变换和邻接矩阵乘法拆开了方便调试。GCNEncoder堆了两层第一层ReLU第二层不加激活。如果你的节点特征维度很高比如超过100可以在两层之间加Dropout比例0.3到0.5防止过拟合。adj_norm需要提前算好不要在forward里重复计算否则训练速度会慢很多。2.3 LSTM时序建模与编码器-解码器结构GCN提取完空间特征后输出的是每个时间步的节点嵌入。接下来要把这些嵌入按时间顺序喂给LSTM。原文用的是编码器-解码器结构编码器部分多个并行GCN模块捕获空间关联然后把带空间相关性的时间序列传入LSTM做时序特征提取编码器生成一个编码向量传给解码器最后全连接层输出预测值。这里有个工程实现上的关键点GCN和LSTM怎么衔接。常见做法有两种一种是先把每个时间步的图快照过GCN得到(N, hidden_dim)的节点嵌入然后把N个节点的嵌入按时间堆叠成(N, T, hidden_dim)再reshape成(N, T, hidden_dim)输入LSTM另一种是把GCN输出直接作为LSTM的初始隐状态。第一种更直观也更容易调试。class GCNLSTM(nn.Module): def __init__(self, num_nodes, in_dim, gcn_hidden, lstm_hidden, pred_len): super().__init__() self.num_nodes num_nodes self.gcn GCNEncoder(in_dim, gcn_hidden, gcn_hidden) self.lstm nn.LSTM( input_sizegcn_hidden, hidden_sizelstm_hidden, num_layers2, batch_firstTrue, dropout0.3 ) self.fc nn.Linear(lstm_hidden, pred_len) def forward(self, x_seq, adj_norm): # x_seq: (N, T, in_dim) N口井T个时间步每个时间步in_dim个特征 N, T, _ x_seq.shape gcn_outs [] for t in range(T): # 每个时间步单独过GCN h self.gcn(x_seq[:, t, :], adj_norm) # (N, gcn_hidden) gcn_outs.append(h) # 堆叠成 (N, T, gcn_hidden) gcn_seq torch.stack(gcn_outs, dim1) # LSTM 时序建模 lstm_out, _ self.lstm(gcn_seq) # (N, T, lstm_hidden) # 取最后一个时间步的输出做预测 last_out lstm_out[:, -1, :] # (N, lstm_hidden) pred self.fc(last_out) # (N, pred_len) return pred这段代码里GCNLSTM的forward先把每个时间步的节点特征过GCN得到空间嵌入序列再整体喂给LSTM。num_layers2表示两层LSTMdropout0.3在层间做正则化。fc层把LSTM最后时间步的隐状态映射到预测长度。如果你的预测任务是输出未来7天水位pred_len7如果只预测下一天pred_len1。注意LSTM的batch_firstTrue输入维度是(batch, seq_len, feature)这里batch维度就是节点数N。3. 数据集构建与两种预测情景的工程实现3.1 时间特征与空间特征的预处理原文用的数据集是56口井、2019年1月1日到2023年12月31日共1826天的观测记录。时间特征包括水温、降雨量、水位空间特征包括经度、纬度、孔口标高、井深。这里有个数据对齐的坑水位和水温是一天记录一次降雨量是一天记录八次。原文做了降采样处理把降雨量按天聚合。常见做法是取日均值或者日累计值具体选哪个要看降雨对地下水位的物理影响机制——如果是入渗补给日累计更合理如果是蒸发影响日均可能更合适。空间特征不随时间变化所以数据集是56×4的静态矩阵。时间特征有两种规模仅水位是1826×56三种时间特征堆叠后是1826×(56×3)。训练集测试集按80/20划分前1460天训练后366天测试。这个划分是时序划分不能随机打乱否则会用未来数据预测过去造成数据泄露。import pandas as pd import numpy as np def preprocess_rainfall(raw_rainfall): raw_rainfall: DataFrame, 每天8条记录 返回日累计降雨量 daily raw_rainfall.resample(D).sum() return daily def build_time_features(water_level, water_temp, rainfall): 三个DataFrameindex为日期columns为井ID 返回 (T, N, 3) 的numpy数组 # 对齐索引 common_idx water_level.index.intersection(water_temp.index).intersection(rainfall.index) wl water_level.loc[common_idx].values # (T, N) wt water_temp.loc[common_idx].values rf rainfall.loc[common_idx].values # 堆叠成 (T, N, 3) features np.stack([wl, wt, rf], axis-1) return features def train_test_split_time(features, train_ratio0.8): 时序划分不能打乱 T features.shape[0] split int(T * train_ratio) return features[:split], features[split:]preprocess_rainfall用resample(D).sum()做日累计如果你的数据是小时级改成resample(D).sum()就行。build_time_features先把三个特征的索引对齐再堆叠成(T, N, 3)。train_test_split_time按时间顺序切分前80%训练后20%测试。注意这里没有做归一化实际训练前要对每个特征做Z-score标准化否则LSTM的梯度会不稳定。3.2 仅空间水位情景适合数据稀缺场景原文把预测情景分成两种仅考虑空间特征和水位特征以及考虑空间特征和三种时间特征。第一种情景适合数据记录不完整的中小型城市只有水位数据也能跑。网络结构上GCN部分不变LSTM的输入维度从3降到1。这种情景的优点是数据要求低很多老旧观测井只有水位记录没有水温传感器更别说降雨量了。缺点是忽略了季节性和气候因素的影响在季节性变化明显的地区预测精度会打折扣。原文的实验结果也印证了这一点仅空间水位情景下MAE和RMSE普遍比三种时间特征情景高0.1左右R²低0.1左右。# 仅空间水位情景的模型配置 model_spatial_only GCNLSTM( num_nodes56, in_dim1, # 只有水位一个特征 gcn_hidden64, lstm_hidden128, pred_len1 ) # 训练循环示例 optimizer torch.optim.Adam(model_spatial_only.parameters(), lr1e-3) criterion nn.MSELoss() for epoch in range(200): model_spatial_only.train() optimizer.zero_grad() pred model_spatial_only(train_x, adj_norm) loss criterion(pred, train_y) loss.backward() optimizer.step() if epoch % 20 0: print(fEpoch {epoch}, Loss: {loss.item():.6f})in_dim1表示只用水位特征。gcn_hidden64和lstm_hidden128是经验值如果你的节点数更多比如超过100可以适当增大。学习率1e-3是Adam的常用起点如果loss不下降降到1e-4试试。训练200轮是个保守值实际要看验证集loss什么时候不再下降早停策略能省不少时间。3.3 空间三种时间特征情景精度优先的配置第二种情景把水温、降雨量、水位三个时间特征都用上LSTM输入维度变成3。原文数据显示这种情景下MAE和RMSE平均比第一种低0.107和0.1049R²高0.1083。代价是数据要求高需要三种传感器都齐全而且数据预处理更复杂。这里有个容易翻车的地方三种特征的量纲差异很大。水位可能是几十米水温是十几度降雨量是几毫米到几百毫米。如果不做标准化直接喂给LSTM降雨量的数值波动会主导梯度更新水位特征被淹没。常见做法是对每个特征单独做Z-score标准化即减去均值除以标准差。注意均值和标准差只能用训练集算然后应用到测试集否则会引入未来信息。def normalize_features(train, test): train, test: (T, N, C) 对每个通道单独做Z-score mean train.mean(axis(0, 1), keepdimsTrue) # (1, 1, C) std train.std(axis(0, 1), keepdimsTrue) 1e-8 train_norm (train - mean) / std test_norm (test - mean) / std return train_norm, test_norm, mean, std # 使用示例 train_norm, test_norm, mean, std normalize_features(train_features, test_features) model_full GCNLSTM( num_nodes56, in_dim3, # 水位、水温、降雨量 gcn_hidden64, lstm_hidden128, pred_len1 )normalize_features里的axis(0, 1)表示在时间维和节点维上算均值和标准差保留通道维。keepdimsTrue保证广播时维度对齐。1e-8是防止除零。标准化后的数据再输入模型训练会稳定很多。如果你发现某个特征的方差特别小比如水温常年波动不超过2度标准化后数值会被放大这时候可以考虑对该特征不做标准化或者用Min-Max归一化替代。4. 训练避坑与调参那些文档没写的血泪经验4.1 邻接矩阵稀疏化别让全连接图拖垮训练56口井如果全连接邻接矩阵有56×563136个元素其中大部分是非零的。GCN每层都要做矩阵乘法计算量随节点数平方增长。如果你的研究区域有几百口井全连接图会让训练慢到无法接受。常见做法是稀疏化只保留每个节点的Top-K个最相似邻居其余置零。K一般取5到10具体看井的空间分布密度。def sparsify_adjacency(A, top_k8): 只保留每个节点的top_k个最大权重邻居 N A.shape[0] A_sparse np.zeros_like(A) for i in range(N): # 排除自身 row A[i].copy() row[i] 0 # 取top_k索引 top_indices np.argsort(row)[-top_k:] A_sparse[i, top_indices] row[top_indices] # 对称化 A_sparse (A_sparse A_sparse.T) / 2 return A_sparsesparsify_adjacency对每一行取Top-K然后对称化保证无向图。top_k8是个经验值井分布密集的区域可以取小一点稀疏区域取大一点。稀疏化后记得重新归一化否则GCN的传播会不稳定。4.2 损失函数选择MSE之外还有哪些选项地下水位预测是回归任务MSE是最常用的损失函数。但MSE对异常值敏感如果某口井的水位记录有跳变传感器故障或人工记录错误MSE会被拉偏。这时候可以考虑Huber损失它在误差较小时等价于MSE误差较大时等价于MAE对异常值更鲁棒。# Huber损失 criterion nn.HuberLoss(delta1.0) # 或者自定义加权MSE对枯水期样本加大权重 class WeightedMSELoss(nn.Module): def __init__(self, weight_low2.0): super().__init__() self.weight_low weight_low def forward(self, pred, target, is_low_season): loss (pred - target) ** 2 weight torch.where(is_low_season, self.weight_low, 1.0) return (loss * weight).mean()HuberLoss的delta控制MSE和MAE的切换点默认1.0。WeightedMSELoss是个自定义示例对枯水期样本加大权重因为枯水期水位变化对水资源管理更关键。is_low_season是个布尔张量标记每个样本是否属于枯水期。这个权重可以根据业务需求调整丰水期权重1.0枯水期2.0到3.0。4.3 常见问题排查现象一训练loss震荡剧烈不收敛。原因通常是学习率太大或者邻接矩阵归一化方式不对。解决把学习率从1e-3降到1e-4检查邻接矩阵是否做了对称归一化确保每行和为1。现象二验证集loss远高于训练集loss过拟合。原因可能是模型参数太多或者训练数据太少。解决增大Dropout比例到0.5减少GCN隐藏层维度或者加L2正则化。如果数据量确实少考虑用数据增强比如对水位序列加小噪声。现象三某些井的预测误差特别大其他井正常。原因可能是这些井的空间位置孤立邻接矩阵里它们的邻居权重很低GCN提取不到有效空间特征。解决检查这些井的坐标如果确实远离其他井考虑在邻接矩阵里给它们加自环权重或者单独用单井模型处理。现象四预测值滞后于真实值波峰波谷对不上。原因通常是LSTM的序列长度不够或者GCN和LSTM的衔接方式有问题。解决增大输入序列长度从7天增加到14天或30天检查GCN输出是否按时间顺序正确堆叠别把时间步搞混了。现象五训练时GPU显存溢出。原因可能是邻接矩阵太大或者batch size太大。解决稀疏化邻接矩阵减小batch size或者用梯度累积模拟大batch。如果节点数超过500考虑用图采样方法每个batch只采样部分节点和边。5. 从复现到落地模型验证与工程化技巧5.1 评价指标解读MAE、RMSE、R²怎么用原文用了三个指标MAE、RMSE、R²。MAE是平均绝对误差单位和水位一样直观反映预测偏差的平均水平。RMSE是均方根误差对大误差更敏感如果RMSE远大于MAE说明存在个别大偏差样本。R²是拟合系数越接近1越好但要注意R²对异常值也敏感而且如果测试集的水位方差很小R²会偏低这不代表模型差。实际工程中我一般会额外看两个指标最大绝对误差MaxAE和预测区间覆盖率PICP。MaxAE告诉你最坏情况下模型偏多少PICP告诉你预测区间是否可靠。如果做防灾减灾MaxAE比MAE更重要因为极端情况下的偏差可能导致误判。def evaluate(y_true, y_pred): y_true, y_pred: (N, T) mae np.mean(np.abs(y_true - y_pred)) rmse np.sqrt(np.mean((y_true - y_pred) ** 2)) ss_res np.sum((y_true - y_pred) ** 2) ss_tot np.sum((y_true - np.mean(y_true)) ** 2) r2 1 - ss_res / (ss_tot 1e-8) max_ae np.max(np.abs(y_true - y_pred)) return {MAE: mae, RMSE: rmse, R2: r2, MaxAE: max_ae}evaluate函数一次性算出四个指标。ss_tot加1e-8防止除零。如果你的数据里有些井的水位几乎不变ss_tot会很小R²可能为负这时候R²参考价值有限重点看MAE和MaxAE。5.2 多井同步预测的输出组织GCN-LSTM的一个优势是一次性输出所有井的预测值不需要为每口井单独跑模型。但输出组织上有个细节模型输出的是(N, pred_len)的张量N是井数pred_len是预测天数。如果你要按井查看结果需要转置成(pred_len, N)再按井ID索引。def organize_predictions(pred_tensor, well_ids): pred_tensor: (N, pred_len) well_ids: list of well identifiers 返回 DataFrame, index为日期, columns为井ID pred_np pred_tensor.detach().cpu().numpy() # (N, pred_len) pred_df pd.DataFrame(pred_np.T, columnswell_ids) return pred_dforganize_predictions把模型输出转成DataFrame方便后续分析和可视化。pred_np.T把(N, pred_len)转成(pred_len, N)列名是井ID。如果你要对比真实值和预测值把真实值也转成同样格式直接做差就行。5.3 模型保存与加载的注意事项PyTorch保存模型有两种方式保存整个模型torch.save(model, path)和只保存参数torch.save(model.state_dict(), path)。推荐后者因为保存整个模型会依赖具体的类定义换环境容易加载失败。加载时先实例化模型结构再load_state_dict。# 保存 torch.save(model.state_dict(), gcn_lstm_weights.pth) # 加载 model_loaded GCNLSTM(num_nodes56, in_dim3, gcn_hidden64, lstm_hidden128, pred_len1) model_loaded.load_state_dict(torch.load(gcn_lstm_weights.pth)) model_loaded.eval()保存时只存state_dict加载时先建模型再加载参数。model.eval()切换到推理模式关闭Dropout和BatchNorm的更新。如果你在GPU上训练、CPU上推理加载时加map_locationcpu。5.4 一个容易忽略的细节邻接矩阵的归一化时机邻接矩阵的归一化要在训练前做好不要在forward里重复算。而且归一化后的矩阵要保存下来训练、验证、测试用同一个。如果每次forward都重新归一化计算量大不说还可能因为数值精度问题导致训练不稳定。def normalize_adjacency(A): 对称归一化: D^{-1/2} A D^{-1/2} N A.shape[0] D np.diag(A.sum(axis1)) D_inv_sqrt np.linalg.inv(np.sqrt(D 1e-8)) A_norm D_inv_sqrt A D_inv_sqrt return A_norm # 训练前算一次 A_norm normalize_adjacency(A_sparse) A_norm_tensor torch.FloatTensor(A_norm)normalize_adjacency做对称归一化1e-8防止除零。算完后转成Tensor训练时直接传入。如果你的图是无向的A应该是对称矩阵归一化后仍然对称。如果A不对称检查构建邻接矩阵时是否做了对称化处理。从那以后我每次复现GCN-LSTM类模型都强制走一遍「邻接矩阵检查 → 特征标准化 → 时序划分验证 → 损失曲线监控」这四步少一步都可能在后头翻车。希望帮到你。本文还有配套的精品资源点击获取