简介基于GCN-LSTM的区域级多井地下水位时空预测模型PDF文档面向环境科学专业学生、水务工程技术人员及地下水研究相关领域人员用于解决多井地下水位同步预测与时空特征提取问题。文档以成都56口国有观测井五年数据为案例系统讲解观测井图结构构建、空间自相似矩阵与属性自相似矩阵计算、GCN与LSTM组合网络结构、编码-解码器搭建及实验对比对降雨量降采样、训练集测试集划分等细节也有说明便于复现和扩展应用。资源包共1个PDF文件大小约3.47MB内容完整聚焦模型原理与实现目前已有419人学习。通过这份文档可掌握GCN-LSTM在区域级地下水预测中的建模思路也能为水务管理、防灾减灾及复杂动态环境下的地下水流动研究提供参考。1. GCN-LSTM 时空预测模型为什么单井模型会漏掉空间关联做地下水位预测的人通常第一反应是上 LSTM 或者 Transformer 直接跑单井历史序列但换成区域级多井场景后很容易翻车你把 56 口井当成 56 个独立任务去拟合训练出来的模型在井与井之间毫无信息交互遇到观测井数据缺失或者水位联动明显的区域预测曲线就跟不上真实波动的拐点。GCN-LSTM 的做法完全不同——先把观测井按空间位置和属性相似度组织成一张图用图卷积网络挖井间空间关联再用 LSTM 挖时间关联最后一次性输出全部井的水位预测值。这份资源来自一篇硕士论文数据是成都市 56 口国有观测井 2019—2023 五年共 1826 天的真实记录模型结构、数据预处理流程和完整的性能对比表都在里面适合环境科学专业做毕业设计的学生、水务工程的技术人员以及想把图神经网络落到水文时序预测上的研究者。2. 模型原理与图结构构建从邻接矩阵到两层 GCN 的前向传播2.1 为什么必须建图井与井之间的水位本来就是连通的地下水位不是孤立变量。同一区域内的观测井井间距可能只有几公里含水层是连通的一口井的抽水、补给会直接影响邻井的水位。如果建模时把每口井当作独立样本等于人为切断了这种物理上的连通性。GCN 解决的就是这个问题把 N 口井看作图上的 N 个节点井与井之间的边表示空间关联强度节点属性是井的经度、纬度、孔口标高、井深这四个不随时间变化的空间特征再加上随时间变化的水位、水温、降雨量序列就构成了一张完整的时空图。这里有个关键设计图中节点之间的边权重不是随便给的而是由空间权重矩阵 W 和属性自相似矩阵 C 相乘得到。W 描述两口井在图结构中的空间位置相似性C 描述单口井与其它井的属性相似情况二者逐元素相乘后再参与图卷积运算相当于把「离得近」和「长得像」两个维度的信息同时编码进邻接关系里。2.2 空间权重矩阵 W距离平方反比公式及其物理含义空间权重矩阵的计算公式是论文里给出的核心公式之一import numpy as np def build_spatial_weight_matrix(coords): 构建空间权重矩阵 W coords: array of shape (N, 2)每行是 (经度, 纬度) 返回 W: array of shape (N, N) n len(coords) W np.zeros((n, n)) for i in range(n): for j in range(n): if i ! j: # 欧式空间距离假设井间距离远大于海拔差忽略高程影响 d np.linalg.norm(coords[i] - coords[j]) W[i, j] 1.0 / (d ** 2) return W这段代码的逻辑是井 i 和井 j 之间的空间距离越近权重越大对水位的影响越强。取距离平方的反比而不是直接取反比是为了突出近距离井之间的强关联——距离翻一倍权重降到原来的四分之一。论文里特别说明了一个近似井与井之间的距离 d 远大于海拔差 h所以两口井连线与水平面的夹角可以忽略不计。实际处理时如果数据里井间距确实很小而海拔差很大这个假设就不成立了需要改用三维距离这是后话。2.3 属性自相似矩阵 C四个空间特征做余弦相似度空间权重矩阵只考虑了地理位置但两口井即使离得很近如果井深和孔口标高差异大水位响应特征也会不同。属性自相似矩阵就是用来补这个维度的。每口井有四个空间特征向量 v (经度, 纬度, 孔口标高, 井深)任意两口井之间的属性相似度用余弦相似度计算def build_attribute_similarity_matrix(features): 构建属性自相似矩阵 C features: array of shape (N, 4)每行是 (经度, 纬度, 孔口标高, 井深) 返回 C: array of shape (N, N)C[i,j] cos(theta) n len(features) C np.zeros((n, n)) for i in range(n): for j in range(n): v1 features[i] v2 features[j] cos_theta np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2)) C[i, j] cos_theta return C注意这里有个细节余弦相似度对向量模长不敏感只关心方向。这意味着孔口标高 500 米的井和标高 100 米的井只要四个特征的比例关系接近相似度就高。实际使用时建议先对特征做标准化否则经度纬度的数值范围会压过井深和标高导致余弦值被经纬度主导。最终邻接矩阵 A W * C也就是两个矩阵逐元素相乘。这样构造出来的 A 既包含「井间距」的物理信息又包含「井属性」的静态信息作为 GCN 的输入。论文里还提到对 A 做了归一化处理即对称归一化形式 D⁻¹/² A D⁻¹/²目的是防止节点度数差异导致数值稳定性问题。2.4 GCN 前向传播公式拆解ReLU 到 Softmax 的两层结构论文给出的两层 GCN 前向传播公式是Z softmax(A_hat relu(A_hat X W0) W1)其中 A_hat 是归一化后的邻接矩阵X 是节点特征矩阵N×FW0 是第一层权重F×HW1 是第二层权重H×CC 是分类数或输出维度。第一层 GCN 用 ReLU 激活第二层用 Softmax。在分类任务里Z 的每一行是该节点属于各个类别的概率在回归任务里Softmax 会被替换成线性输出或者套一个全连接层把特征维度映射到 1。import torch import torch.nn as nn class TwoLayerGCN(nn.Module): def __init__(self, in_features, hidden_dim, out_dim, dropout0.5): super().__init__() self.gc1 GraphConv(in_features, hidden_dim) # 第一层: ReLU self.gc2 GraphConv(hidden_dim, out_dim) # 第二层: Softmax/线性 self.dropout nn.Dropout(dropout) def forward(self, x, adj): h self.gc1(x, adj) h torch.relu(h) h self.dropout(h) z self.gc2(h, adj) return torch.softmax(z, dim-1) # 分类任务; 回归任务改为线性输出参数说明in_features 是节点特征维度在论文的场景里是 56 口井×历史水位序列长度或空间特征维度hidden_dim 是隐藏层维度一般取 16 或 32out_dim 是输出维度多井预测时等于井数量 56。损失函数用的是交叉熵对所有带标签的节点计算损失后反向传播更新参数。对这篇论文的场景要特别注意GCN 不是用来做节点分类的而是用来输出 w 每个观测井在下一时刻的预测水位值所以第二层的 Softmax 在实际的预测代码里通常会被去掉换成全连接层或直接输出回归值论文第 5.3.2 节提到的「编码器-解码器结构」才是最终落地形式。3. 数据组织与模型搭建成都 56 口井五年数据的两种输入规模3.1 数据集到底长什么样三个时间特征加四个空间特征论文用的数据集是成都市 56 口国有观测井时间跨度 2019 年 1 月 1 日到 2023 年 12 月 31 日共 1826 天。时间特征有三种水位、水温、降雨量。水位和水温每天记录一次降雨量一天记录八次为了统一维度做了降采样处理——对当天的八次降雨量数据做聚合把一天内的降雨量压缩成一个值。空间特征有四个经度、纬度、孔口标高、井深这些是不随时间变化的静态数据。数据规模分两种空间特征矩阵56×4共 56 口井每口井 4 个空间特征。水位时间特征数据集1826×56共 1826 天、56 口井每口井每天一个水位值。三时间特征数据集1826×(56×3)每口井每天有水位、水温、降雨量三个值总列数为 56×3168 列。这里容易混淆的是矩阵的行列含义。很多初学者会把 1826×56 理解成 1826 个样本、56 个特征但实际上这个矩阵的每一行是一天每一列是一口井在使用时需要做转置或 reshape把井 ID 变成样本维度。在 GCN 里节点是 56 口井每个节点的特征向量是从时间序列里滑窗取出的历史窗口我一般会先按井分组再用滑窗构造样本避免直接把原始 1826×56 矩阵丢进模型。3.2 训练集和测试集怎么切80/20 比例与时间序列的切分陷阱论文明确写了 80% 训练集、20% 测试集即前 1460 天用于训练后 366 天用于测试。这是一个严格按时间顺序的切分不能随机打乱——如果随机切分模型会「偷看」未来数据导致预测指标虚高这种错误在时序预测里非常常见。import pandas as pd from sklearn.preprocessing import StandardScaler def load_well_data(water_path, rain_temp_path, spatial_path): 加载水位、降雨/水温、空间特征数据按时间顺序切分 返回训练集、测试集和标准化器 water pd.read_csv(water_path, index_col0) # shape: 1826 x 56 rain_temp pd.read_csv(rain_temp_path, index_col0) # shape: 1826 x (56*3) # 时间顺序切分: 前80%训练, 后20%测试 split_idx int(len(water) * 0.8) # 1460天 water_train, water_test water.iloc[:split_idx], water.iloc[split_idx:] rain_temp_train, rain_temp_test rain_temp.iloc[:split_idx], rain_temp.iloc[split_idx:] # 空间特征不随时间变化, 不参与切分 spatial pd.read_csv(spatial_path, index_col0) # shape: 56 x 4 return (water_train, water_test), (rain_temp_train, rain_temp_test), spatial参数说明split_idx 就是 1460用 int(len(water) * 0.8) 是为了让切分边界对数据长度变化自动适配。StandardScaler 的 fit 必须只用训练集数据测试集用同一套均值和方差转换不能重新 fit。3.3 GCN-LSTM 网络结构编码器-解码器与两种预测场景模型主体是编码器-解码器结构。编码器部分先把考虑了时空特征的邻接矩阵56×56输入多个并行的 GCN 模块GCN 捕获节点间的空间结构和属性关系挖掘不同观测井之间的空间关联性。GCN 的输出是带有空间相关性的时间序列数据再传给 LSTMLSTM 做序列特征分析和进一步特征提取挖掘时间关联关系。编码器最终生成一个编码向量传给解码器解码器通过全连接层把时空特征转换回原始空间输出每口井的地下水位预测值。这里要理解 GCN 和 LSTM 的衔接方式。GCN 的输入是节点特征矩阵在时间维度上是一个时间片的快照LSTM 处理的是序列所以中间必然有一个 reshape 操作把 GCN 输出的特征按时间步展开组成 (batch, seq_len, hidden_dim) 形状的序列数据传入 LSTM。我见过不少实现直接把 GCN 输出的二维矩阵塞进 LSTM结果维度对不上报错后就开始怀疑模型结构有问题。import torch import torch.nn as nn class GCNLSTMEncoder(nn.Module): def __init__(self, gcn_in_dim, gcn_hidden, lstm_hidden, num_wells): super().__init__() self.gcn TwoLayerGCN(gcn_in_dim, gcn_hidden, gcn_hidden) self.lstm nn.LSTM(input_sizegcn_hidden, hidden_sizelstm_hidden, batch_firstTrue) self.fc nn.Linear(lstm_hidden, num_wells) def forward(self, x_seq, adj): # x_seq: (batch, seq_len, num_wells, feat_dim) batch, seq_len x_seq.shape[0], x_seq.shape[1] gcn_out [] for t in range(seq_len): xt x_seq[:, t, :, :] # (batch, num_wells, feat_dim) h self.gcn(xt, adj) # 对每个时间片做图卷积 gcn_out.append(h) gcn_seq torch.stack(gcn_out, dim1) # (batch, seq_len, num_wells, hidden) # 把井维度合并进 batch, LSTM 对每口井独立编码时间特征 b, s, n, h gcn_seq.shape gcn_seq gcn_seq.permute(0, 2, 1, 3).reshape(b * n, s, h) lstm_out, _ self.lstm(gcn_seq) # (b*n, seq_len, lstm_hidden) last_hidden lstm_out[:, -1, :] # 取最后一个时间步 pred self.fc(last_hidden) # (b*n, num_wells) 不经过 softmax return pred.view(b, n)模型里的几个关键参数gcn_hidden 一般取 16~32太小特征提取能力不足太大容易过拟合lstm_hidden 取 32 或 64要和 gcn_hidden 保持在同一个量级num_wells 就是 56。forward 里最核心的是维度转换GCN 逐时间片处理LSTM 再按井独立编码时序最后通过全连接层输出 56 个节点的水位预测值。论文设计了两种工程场景第一种只输入空间特征和水位特征适用于数据记录不完整的中小型城市第二种输入空间特征加水位、水温、降雨量三个时间特征适用于数据充沛的大型城市。从论文表 5-1 的数据看加入水温降雨量后绝大多数井的 MAE 平均下降 0.107、RMSE 平均下降 0.1049、R² 平均提升 0.1083效果提升明显但数据采集成本和预处理复杂度也上去了。4. 从复现到迁移指标评估、调参经验与自备数据接入4.1 三个评价指标怎么解读MAE、RMSE、R² 的合理取值范围论文用了三个回归指标MAE平均绝对误差、RMSE均方根误差、R²决定系数。MAE 衡量预测值和真实值之间的平均绝对偏差单位和水位一样是米RMSE 对大的误差更敏感因为先平方再开方会把个别离谱的预测点放大R² 反映模型对真实值方差的解释程度越接近 1 越好。从论文表 5-1 看56 口井在「空间水位水温降雨量」场景下R² 基本分布在 0.83~0.97 之间表现最好的 well-13 达到 0.9634MAE 最低的 well-49 只有 0.0567。需要注意R² 为 0.97 不代表误差只有 3%它只说明模型解释了 97% 的方差实际误差要看 MAE 和 RMSE 的具体值。对地下水位预测来说MAE 在 0.1 米以内已经算相当不错0.3 米以上的预测结果就需要警惕是不是数据质量或模型结构出了问题。4.2 训练参数怎么调学习率、批次大小、序列长度的经验值论文没有展开训练超参数但按 GCN-LSTM 时空预测的常见做法我一般会这么设学习率 0.001 起步Adam 优化器训练 200~300 轮后 loss 不降就降一半学习率。GCN 收敛比纯 LSTM 慢学习率大了容易震荡小了收敛太慢。批次大小 32 或 64。56 口井本身节点数不多批次主要作用在时间窗口维度batch 太大容易让模型记住训练集的均值。序列长度时间窗口取 7 天或 30 天。论文的数据是逐日采样7 天窗口捕捉一周内的水位波动规律30 天窗口能覆盖更长的退水过程。窗口太长会导致训练样本数锐减因为总天数只有 1826 天。Dropout 设置在 0.3~0.5 之间。GCN 层和 LSTM 层之间各加一个 dropout防止模型过度依赖某几口井的特征。4.3 迁移到自己的数据把成都 56 井换成任意井群数据拿到这份资源后很多人会用在自己的区域数据上。替换数据的步骤很清楚把成都 56 口井的空间特征 CSV 换成你的井群坐标和井深信息用 2.2 和 2.3 里的代码重新计算 W 矩阵和 C 矩阵时间特征 CSV 换成你的水位、水温、降雨量序列。这里有几个容易忽略的点井数量变了邻接矩阵维度要跟着变。假设你的区域只有 20 口井W 和 C 都是 20×20模型输出的 num_wells 也要改成 20不能直接用 56 口井训练好的权重。经纬度坐标要统一投影。如果有的井是经纬度、有的是平面坐标直接算欧氏距离会得到荒谬的结果建议统一转成 UTM 投影后再算距离。属性相似度矩阵的特征要事先归一化。论文里直接用了原始值做余弦相似度但不同地区的井深范围差异很大比如成都的井深可能是 50~200 米换到北方平原可能变成 200~500 米不归一化会导致 C 矩阵被井深这个维度主导。5. 避坑指南GCN-LSTM 复现里的五个典型翻车点5.1 邻接矩阵归一化方式用错梯度消失或数值爆炸现象训练 loss 一开始就变成 NaN或者模型输出全是同一个值。原因直接用了原始 W*C 矩阵作为 GCN 的邻接矩阵没有做归一化。原始矩阵每一行的权重和差异很大节点度数高的井在聚合邻居特征时数值会远大于其它井经过两层图卷积后数值爆炸。解决改用对称归一化 D⁻¹ᐟ² A D⁻¹ᐟ²也就是先算度矩阵 D对角矩阵每个对角元素是 A 对应行的和再对每个元素按公式变换。我一般直接用 PyTorch Geometric 里的 GCNConv 层它内部自带归一化省去手动实现的出错风险。5.2 降雨量降采样方式错误把 8 次观测直接删掉 7 次现象加入降雨量特征后预测效果反而变差R² 明显下降。原因降雨量一天记录八次为了维度统一直接取第一个值或最后一个值丢失了降雨过程的强度信息。一天的降雨分布不均匀一次短时强降雨可能集中在某一个小时取首值或末值都会引入噪声。解决论文说做了降采样处理但没有明确说是取均值还是求和。我建议取每天的累计降雨量8 次求和因为累计值保留了「这天下了多少雨」的核心信息而取均值会把降雨强度平均稀释掉。注意如果原始记录里单位是毫米/次求和后单位还是毫米不需要额外转换。5.3 GCN 输出维度和 LSTM 输入维度接不上现象运行到 GCN 输出转 LSTM 输入时报错比如 Expected 3D tensor but got 2D tensor。原因GCN 输出的形状是 (batch, num_wells, hidden_dim)但 LSTM 要求输入是 (batch, seq_len, input_size)中间缺少时间维度的展开。解决在 GCN 和 LSTM 之间加一个维度重排的代码把每个时间步的 GCN 输出按时间顺序堆叠成序列再把井维度合并进 batch 维度。第 3.3 节代码里已经给出了这种处理方式核心是 gcn_seq.permute(0, 2, 1, 3).reshape(b * n, s, h)这个 reshape 操作很多人第一次写都会漏。5.4 用全量数据训练和论文结果对不上现象复现出来的 MAE 比论文表 5-1 里的大很多或者 R² 低 0.1 以上。原因很可能把 1826 天全部数据拿去训练了测试时没有按照 80/20 的时间顺序切分。时序预测和普通回归不同样本之间有强时间依赖如果训练集里混入了测试时间段的数据模型在训练时就见过了未来的水位测试指标会异常好但如果反过来把测试集放在训练集前面模型又完全没见过未来模式指标会异常差。解决严格用前 1460 天训练、后 366 天测试中间不能有任何随机打乱。标准化时也只对训练集拟合均值和方差测试集直接用训练集的参数做变换。5.5 空间特征参与标准化后被改变物理含义现象属性相似度矩阵 C 算出来的值全部接近 1模型无法区分不同井的差异。原因经度纬度的数值范围在 100 左右孔口标高可能在几百米井深可能在几十米四个特征不在一个量级。原始值直接算余弦相似度时经度和纬度因为数值大而主导了相似度井深和标高的影响被淹没。解决先在空间特征列上做标准化但注意标准化后特征的物理含义就变成了「偏离均值的程度」不再是原始值。如果想让 C 矩阵保持可解释性可以按特征分组做缩放经纬度归一到 0~1标高和井深各自归一到 0~1再计算余弦相似度。6. 验证模型效果的正确姿势从 well-56 对比到波形拐点检查拿到训练好的 GCN-LSTM 模型第一件事不是看总指标而是挑几口有代表性的井画「预测值 vs 真实值」的曲线。论文选了 well-1、well-10、well-21、well-31、well-38、well-53 这六口井做了可视化展示还特别对比了 well-56 在 GCN-LSTM 和单井 STA-LSTM 两种模型下的表现。这个思路值得直接用在自己的数据上well-56 在 GCN-LSTM 下的 MAE 是 0.02415STA-LSTM 是 0.02437RMSE 分别是 0.02568 和 0.02857R² 分别是 0.97472 和 0.96673。指标差距虽小但结合曲线看GCN-LSTM 在波峰和波谷的捕捉上明显更准转折点没有滞后。具体验证分三步。第一步取测试集最后 50 天的数据画出预测值和真实值的对比折线重点观察曲线的峰值和谷值是否对齐。如果整体趋势一致但峰值偏低说明模型对极端水位变化不敏感可以考虑增加 GCN 层数或扩大 LSTM 隐藏维度。第二步计算每口井的逐日误差分布找出误差最大的几口井检查它们的空间特征是否有异常比如井深特别小或靠近河流的井水位波动可能被局部水文条件主导而不是邻井关联主导。第三步用「加入水温降雨量 vs 不加」两个版本分别预测同一批数据如果加了更多特征的版本在某些井上反而更差说明这部分特征在局部区域是噪声需要特征筛选。最后说说我的固执习惯每次跑完一个版本的模型我都会强制保存三样东西——训练集和测试集的切分索引、标准化器的均值和方差、最后五个 epoch 的 loss 曲线。原因很简单切分索引决定了实验结果能否复现标准化参数决定了模型能否部署到新数据上loss 曲线决定了你是欠拟合还是压根没收敛。这三样缺一样前面花了几天调的参数等于白调。从那以后我每次复现类似的水文时空预测模型都会强制走一遍这个流程。希望帮到你。本文还有配套的精品资源点击获取