区域地下水位预测的GCN-LSTM实战：从邻接矩阵构建到PyTorch实现

简介一份聚焦区域级多井地下水位时空预测的学术文档面向环境科学、水务工程及机器学习交叉领域的研究者。内容提出融合图卷积网络与长短期记忆网络的预测模型通过构建基于观测井位置与属性的空间邻接矩阵将成都五十六口国有井五年日观测数据组织成图结构结合水位、水温、降雨时序特征实现多井水位同步预测。文中详细给出图卷积传播公式、空间自相似矩阵与属性自相似矩阵的构造方法并对比仅用空间特征和水位特征与加入温度、降雨等气象因素两种情景的预测效果说明引入气象信息后可提升预测精度。资源为单个PDF文件体积3.47MB内容包含完整模型原理、网络结构图和实验分析便于直接阅读与复现。已有419人学习浏览适合需要借鉴时空图神经网络在水文预测中应用方案的研究人员与技术人员。1. 区域级地下水位预测为什么要上 GCN-LSTM区域水资源管理里最头疼的不是单口井的水位预测而是把全县、全流域几十上百口监测井的水位同时预报出来。传统做法是每口井各自建立 ARIMA、SVR 或者单站 LSTM 模型这等于默认井与井之间没有关系。可实际上抽水降深会在含水层里传播河流补给只影响沿岸井断层和岩性变化会切断水力联系——这些空间上的相互影响单井模型完全没接住。GCN-LSTM 要解决的就是这个问题GCN 把监测井网当成一张图沿着井间连接关系聚合邻居信息刻画空间依赖LSTM 沿着时间轴捕捉每口井的滞后响应刻画时间依赖。二者组合起来就是一个同时带空间和时间结构的预测模型。这篇按可复现的标准走一遍完整流程从数据格式、邻接矩阵构建到模型训练、验证和避坑让读者拿自己的井网数据也能跑出可用结果。2. 井网为什么是一张图GCN 能捕捉到的水力联系2.1 地下水位数据里的空间依赖长什么样地下水位不是一个孤立的数值序列。同一含水层里一口井抽水周围观测井的水位会在几天到几个月后出现不同程度的下降下降幅度和滞后时间取决于井间距、含水层导水系数和储水系数。河流补给带附近的井水位跟着河水位涨落山前补给区的井水位受降雨入渗影响更大。这些关系有一个共同特点它们可以被“井间距离 水文地质条件”近似表达而这种成对关系天然适合用图来建模。如果把每口监测井当作图上的一个节点两井之间有水力联系就连一条边那么一次水位变化就是信息在图上传播的过程。GCN图卷积网络做的事情就是让每个节点在每一层卷积里把邻居节点的特征加权求和再经过线性变换和非线性激活逐层把“邻居的信息”揉进自己的表征里。区域级井网通常有几十到几百个节点节点数不大边也稀疏GCN 在这种规模下计算开销很低比 Transformer 这类需要全局注意力的结构省得多。2.2 为什么是 GCNLSTM 的组合而不是纯 GCN 或纯 LSTM纯 LSTM 方案最容易想到也最容易实现但它有一个原则性问题每口井的 LSTM 参数是独立训练还是共享参数如果独立训练n 口井就是 n 个模型数据量被切碎井多的区域还能用井少的区域很容易过拟合如果共享参数所有井的时间序列被当成同一分布采样但山区井和河谷井的动态规律并不相同模型会被不同分布拉扯。更关键的是纯 LSTM 的输入特征里如果不显式加入邻居井的水位它就完全没有关于空间依赖的信息。纯 GCN 方案把每个时刻的井网水位作为图信号做卷积能捕捉空间同步的相关性但地下水位响应是滞后的——抽水井的水位变化传到观测井需要时间。用纯 GCN 需要把这个滞后效应手工做成多阶邻居特征阶数很难确定物理意义也不直观。GCN-LSTM 的组合思路是把问题拆成两步在每一个时间步上先用 GCN 做一次空间聚合让每口井的当前特征携带邻居井的同步信息再把聚合后的特征序列交给 LSTM让 LSTM 去学习单井的时间滞后规律。空间和时间分开处理每个模块的任务都更单纯参数也更好调。2.3 空域 GCN 的传播规则用 PyTorch 实现一个 GCN 层图卷积有两大流派。谱域 GCN 需要计算拉普拉斯矩阵的特征分解图一变特征分解就得重算区域级井网的条件一变比如增加一口新井整个模型都要重新训练很不实用。空域 GCN 直接定义邻居聚合操作不用做特征分解实现简单、训练稳定是井网场景下更常用的做法。它的核心公式是H ReLU(D⁻¹/² Â D⁻¹/² H W)其中 Â A I 是加了自环的邻接矩阵D 是 Â 的度矩阵W 是可学习的权重矩阵。左侧的 D⁻¹/² Â D⁻¹/² 叫对称归一化邻接矩阵作用是让度数高的节点邻居多的井不会被累积求和放大数值。代码实现很直接import torch import torch.nn as nn import torch.nn.functional as F class GCNLayer(nn.Module): def __init__(self, in_features, out_features): super().__init__() self.linear nn.Linear(in_features, out_features) def forward(self, x, adj_norm): # x: [N, F_in]当前时刻所有井的特征矩阵 # adj_norm: [N, N]已做对称归一化的邻接矩阵 x torch.mm(adj_norm, x) # 邻居特征聚合 return F.relu(self.linear(x)) # 线性变换 激活输入里adj_norm是提前算好、训练过程中不再变化的矩阵x是某一时刻所有井的特征包括水位、降雨、蒸发等变量。torch.mm(adj_norm, x)这一步就是把邻居井的特征按边的权重累加到当前井上然后过一层线性变换。GCN 层数一般取 2 层就够堆太多层会出现过平滑现象——所有节点的表征趋同井间的差异性被磨平这在节点数只有几十个的区域级场景里尤其明显。自环很重要加了自环节点在聚合时才会保留自己上一层的特征。如果不加第一层卷积后每个节点就只能看到邻居信息自己的历史水位信息反而丢了。对称归一化则是为了平衡节点度数的差异否则靠近城市、井网密集的区域数值会被放大。3. 数据准备与邻接矩阵构造可复现的关键一步3.1 字段、时间粒度和目标变量的统一GCN-LSTM 的输入是三维以上的张量比普通表格数据多了一个图结构维度所以数据准备的第一步是把原始数据整理成规整的矩阵形式。以月尺度预测为例需要的核心字段如下字段说明常见单位date观测日期年-月-well_id监测井编号-gw_level地下水位埋深或高程mprecip月累计降雨量mmevap月累计蒸发量mmabstraction月开采量如有10⁴ m³river_level附近河流水位如有m这里第一个容易踩坑的是“埋深”和“高程”两个口径。水位埋深是地面到地下水面的距离数值越大说明水位越低高程是绝对标高数值越大水位越高。两种口径都能做预测但要注意如果目标井里有几口用的是埋深、几口用的是高程混合输入会让模型学到虚假的“井间差异”。统一做法是全部转成同一基准或者至少对每一口井做独立的标准化后再进入模型。时间粒度上区域地下水监测数据以月报为主日尺度数据在欠发达地区很难拿全。用月尺度时lookback 窗口通常取 12 到 24 个月这覆盖了一个完整水文年内的丰枯周期。数据集时长方面至少要有 8 到 10 年的历史数据才值得上深度学习模型数据少于 5 年时优先考虑线性回归或轻量梯度提升把 GCN-LSTM 留到数据量足够的场景。3.2 邻接矩阵三种构造方式距离阈值、KNN 和水文分区图上边的定义直接决定 GCN 能聚合到什么信息。常见的做法有三种各有取舍。第一种是距离阈值法两井之间距离小于某个阈值就建边。优点是无脑、不需要额外资料缺点是阈值敏感——设小了图太稀疏设大了把不同含水层的井也连上了。第二种是 KNN每口井固定找最近的 k 口井建边。优点是每口井至少有一定数量的邻居图结构稳定缺点是会把本来没有水力联系的井硬拉在一起。第三种是结合水文地质分区根据含水层类型、地下水流场、断层位置等先验知识只在同一水文地质单元内的井之间建边。这是最接近物理实际的做法但需要区域水文地质图和相关专家经验资料门槛高。实际项目中我一般先用 KNN 快速跑通取得基线结果后再用手里的水文地质资料把明显的错误边删掉。构建邻接矩阵的代码可以这样写import numpy as np from scipy.spatial.distance import cdist def build_adjacency(coords, methodknn, k5, dist_threshNone): n len(coords) dist cdist(coords, coords) # [N, N] 井间距离矩阵 adj np.zeros((n, n)) if method knn: for i in range(n): idx np.argsort(dist[i])[1:k1] # 最近 k 个邻居跳过自身 adj[i, idx] 1.0 adj[idx, i] 1.0 # 保证对称 elif method distance: if dist_thresh is None: dist_thresh np.percentile(dist, 25) adj[dist dist_thresh] 1.0 elif method gauss: if dist_thresh is None: dist_thresh np.percentile(dist, 25) sigma dist_thresh / 2.0 adj[dist dist_thresh] np.exp(-dist[dist dist_thresh]**2 / (2 * sigma**2)) np.fill_diagonal(adj, 1.0) # 自环 return adj代码里有几个值得注意的点。KNN 方法中argsort返回从小到大排列的索引[1:k1]跳过自身取最近 k 个邻接矩阵必须做对称处理因为 GCN 的归一化公式默认邻接矩阵是对称的否则不对称会产生偏差。距离阈值法里用 25 分位数作为默认阈值这样大约四分之一较近的井对会连边算是一个经验默认值。高斯核方法把硬边变成软权重距离越近权重越大适合水力联系随距离衰减明显的冲积平原地区。有了邻接矩阵还必须做对称归一化否则度数大的节点在聚合时数值会偏大。归一化的实现如下def normalize_adj(adj): deg adj.sum(axis1) deg_inv_sqrt np.power(deg, -0.5) deg_inv_sqrt[np.isinf(deg_inv_sqrt)] 0.0 return deg_inv_sqrt[:, None] * adj * deg_inv_sqrt[None, :]这段代码对应公式 D⁻¹/² Â D⁻¹/²先计算每个节点的度数取 -0.5 次方再左乘右乘。孤立节点度数为 0-0.5 次方会变成无穷大必须替换成 0。归一化后的矩阵在模型训练前转换成torch.float32张量全程不需要参与梯度更新。3.3 时间窗口切分滑动窗口生成样本整理好的数据是一张大表需要切成模型能吃的样本对。输入是过去 lookback 个月的逐月特征输出是未来 horizon 个月的每个井水位。以 lookback12、horizon1 为例def make_samples(features, targets, lookback12, horizon1): # features: [T, N, F]targets: [T, N] T len(features) X, Y [], [] for t in range(lookback, T - horizon 1): X.append(features[t - lookback:t]) Y.append(targets[t horizon - 1]) return np.array(X), np.array(Y)输出X的形状是[样本数, lookback, N, F]Y的形状是[样本数, N]。这里有个常被忽略的细节t从lookback开始因为预测thorizon时刻的水位用的输入只能是t-lookback到t-1的数据绝不能包含当前时刻 t 本身的值。如果监测数据有缺失可以先把缺失月份标记出来构造样本时跳过或者在训练时用掩码把缺失位置从损失函数里剔除。数据集的划分不能用随机切分。地下水位的年周期和多年趋势意味着相邻月份的样本高度相关随机划分会把相互依赖的样本同时放进训练集和测试集造成验证指标虚高。按时间顺序切分是最稳妥的前 80% 的月份做训练后 20% 做验证。如果要做更严格的评估可以用滚动预测的方式每个月重新训练一次模型但计算成本高一般只在最终效果评估时使用。3.4 标准化要在训练段上计算这不是小题大做水位、降雨、蒸发量纲不同数值范围差几个数量级不标准化直接进模型大数值变量会主导梯度。常见做法是 z-score 标准化但必须只在训练段上计算均值和标准差然后用这套统计量去变换验证段和测试段。如果对整个数据集一起算均值方差验证段的分布信息就泄露进了标准化过程这属于一种隐蔽的数据泄露。训练集统计量计算示例train_len int(len(X) * 0.8) x_train, x_val X[:train_len], X[train_len:] mean x_train.mean(axis(0, 1, 2), keepdimsTrue) std x_train.std(axis(0, 1, 2), keepdimsTrue) x_train (x_train - mean) / std x_val (x_val - mean) / stdaxis(0,1,2)是在样本、时间、井三个维度上分别计算每个特征的均值标准差这样每个特征变量各有一套统计量。验证和测试阶段要沿用训练集的 mean 和 std不能重新计算。4. GCN-LSTM 模型实现PyTorch 最小可跑通代码与参数设置4.1 整体结构先空间聚合再时间建模模型的计算顺序是对每个时间步做 GCN 空间聚合得到每个井每个时刻的聚合特征然后把每个井的特征序列送入 LSTM让 LSTM 沿着时间方向更新隐状态最后取 LSTM 最后一个时刻的隐状态经过线性层输出预测水位。这个结构把空间和时间两个维度的计算解耦每一块的输入输出都清晰可查。具体到张量维度输入X是[B, N, T, F]分别对应批大小、井数、时间步数、特征数。GCN 是在每个时间步上独立操作的所以需要把时间维挪到和 Batch 维合并变成[B*T, N, F]做一次图卷积。LSTM 则是每口井独立处理时间序列需要把井维挪到 Batch 维变成[B*N, T, H]。4.2 模型类代码维度变换是唯一的难点import torch import torch.nn as nn class GCNLSTM(nn.Module): def __init__(self, in_channels, hidden_gcn32, hidden_lstm64, num_layers1, out_channels1, dropout0.2): super().__init__() self.gcn GCNLayer(in_channels, hidden_gcn) self.lstm nn.LSTM(hidden_gcn, hidden_lstm, num_layers, batch_firstTrue, dropoutdropout) self.reg nn.Linear(hidden_lstm, out_channels) def forward(self, x, adj_norm): # x: [B, N, T, F] B, N, T, F x.shape x x.permute(0, 2, 1, 3) # [B, T, N, F] x x.reshape(B * T, N, F) # 每个时刻独立做 GCN h self.gcn(x, adj_norm) # [B*T, N, hidden_gcn] h h.reshape(B, T, N, -1) h h.permute(0, 2, 1, 3) # [B, N, T, hidden_gcn] h h.reshape(B * N, T, -1) # 每口井独立过 LSTM out, _ self.lstm(h) # [B*N, T, hidden_lstm] out self.reg(out[:, -1, :]) # 取最后时间步 return out.reshape(B, N, out_channels)这段代码的四个维度变换是整个模型的核心逐一说明x.permute(0, 2, 1, 3)把时间维挪到第二维是为了后面合并 Batch 和时间时各井特征连续排列。reshape(B*T, N, F)合并后每个样本代表“某个时刻的整个井网”GCN 对它做空间卷积。h.reshape(B, T, N, -1)把时间维恢复再permute把井维挪到 Batch 后面reshape(B*N, T, -1)得到的是“某口井的时间特征序列”交给共享参数的 LSTM。out[:, -1, :]取 LSTM 最后一个时间输出的隐状态因为预测目标是未来时刻最后时刻的隐状态包含了整段时间序列的压缩信息。GCNLayer使用的是上一节实现的单层图卷积这里只堆了一层。如果需要两层 GCN可以在__init__里再添加一个GCNLayer并让第二层的输入维度等于第一层输出维度。4.3 训练循环损失函数、优化器和梯度裁剪训练代码用最小化的写法实现便于替换成自己的数据import torch.nn.functional as F from torch.optim import Adam def train_model(model, x_train, y_train, x_val, y_val, adj_norm, epochs200, lr5e-3, devicecuda): model.to(device) opt Adam(model.parameters(), lrlr, weight_decay1e-5) adj_norm adj_norm.to(device) x_tr torch.tensor(x_train, dtypetorch.float32).to(device) y_tr torch.tensor(y_train, dtypetorch.float32).to(device) x_va torch.tensor(x_val, dtypetorch.float32).to(device) y_va torch.tensor(y_val, dtypetorch.float32).to(device) for epoch in range(epochs): model.train() opt.zero_grad() pred model(x_tr, adj_norm) # [B, N, 1] loss F.mse_loss(pred.squeeze(-1), y_tr) loss.backward() nn.utils.clip_grad_norm_(model.parameters(), max_norm1.0) opt.step() if (epoch 1) % 20 0: model.eval() with torch.no_grad(): val_pred model(x_va, adj_norm) val_loss F.mse_loss(val_pred.squeeze(-1), y_va) print(fepoch {epoch1:3d} train {loss.item():.4f} fval {val_loss.item():.4f})参数设置的经验值hidden_gcn32在井数 30 到 100 的区间内够用井数更多时可以加到 64hidden_lstm64对月尺度数据通常足够过大的 LSTM 隐状态在小数据集上容易过拟合lr5e-3配合clip_grad_norm_是稳定的组合直接上1e-2很容易梯度爆炸。LSTM 层数默认取 1除非数据量很大且时间序列有明显的高阶依赖才考虑增加到 2 层。4.4 预测值还原标准化统计量要当成模型的一部分保存模型输出是标准化空间里的预测值不能直接用于业务。还原公式是把 z-score 逆变换def inverse_transform(pred, mean, std): pred pred.squeeze(-1) * std.squeeze(-1) mean.squeeze(-1) return pred这里的mean和std必须是训练集上计算的统计量并且维度要与输出形状对齐。更稳妥的做法是连同标准化统计量、邻接矩阵、井坐标一起保存到文件方便部署阶段直接加载。很多人训练好模型只保存了权重部署时才发现统计量和图结构没存重新算又对不上这是可复现流程里最常见的失误之一。5. 复现避坑指南四个容易让结果翻车的边界问题5.1 数据泄露验证集像是“偷看了”答案现象训练集和验证集的 RMSE 都低得离谱模型表现得比论文里的还好看但换到新数据后完全失灵。原因通常是三种叠加一是数据集划分用了随机shuffle相邻月份样本被分到两边信息互相渗透二是标准化时使用了全体数据的均值标准差验证集的分布特征被带进训练过程三是在构造滑动窗口时输入特征里混入了当前时刻甚至未来时刻的水位。解决方案样本构造严格保证输入是t-lookback到t-1输出是thorizon-1划分按时间顺序切分所有统计量只在训练段上计算最后用时间序列交叉验证比如把 10 年数据分成训练 8 年、验证 2 年复核一遍。这三步都做到数据泄露基本能堵住。5.2 预测出负水位不是模型坏了要先看业务口径现象监测数据里水位埋深全部为正模型预测出一片负数看起来完全不合理。先说一个容易误判的点如果数据用的是水位高程负值是完全合法的华北平原不少地区的地下水高程本来就在零米以下。只有确认业务口径是“埋深”地面往下量恒为正时负数预测才是问题。常见原因有两个。第一某口井的水位序列存在逐年下降趋势z-score 标准化后模型学到的是均值回归遇到趋势性下降时段会在低值区间过冲。第二输出层是线性层没有任何约束训练不充分时中间状态偶尔越界是正常的但成片负值说明模型没有学到合理分布。解决思路如果业务上埋深必然为正给输出层套一个softplus或ReLU让模型输出天然非负更根本的做法是检查数据里是否存在仪器故障、站点迁移导致的跳跃值这种异常会让模型在离群点附近学到错误映射。5.3 井网太稀疏GCN 退化成了每口井自己的线性变换现象加上 GCN 之后模型效果和纯 LSTM 完全一样GCN 部分好像没起作用。先检查邻接矩阵打印每一行的非零个数如果中位数是 1只有自环那 GCN 确实什么都没聚合到。距离阈值法最容易出现这个问题——阈值设得太小只有极近的井对建了边大多数井没有邻居GCN 的聚合操作退化成恒等映射。解决办法是改用 KNN 方法保证每口井至少有 3 到 5 个邻居或者用高斯核软权重让近处井有较大权重远处井有小权重而不是一刀切不连边。另一个常见错误是邻接矩阵构建时忘了对称化导致归一化矩阵不对称GCN 训练不稳定特征被某一方向过度放大。5.4 训练损失震荡、验证指标飘忽小数据集上的不确定性现象同一份代码、同一个随机种子跑了三次结果不一样损失曲线忽高忽低验证 NSE 在 0.3 到 0.7 之间浮动。根本原因是区域级井网数据量太小样本数通常只有一两百个月模型容量稍微大一点就会过拟合。解决策略分级处理先把学习率从5e-3降到1e-3观察损失是否平稳加上梯度裁剪max_norm1.0防止 LSTM 梯度爆炸再不行就把hidden_lstm从 64 减到 32LSTM 层数保持 1 层。最后要固定随机种子保证每次实验可复现def set_seed(seed42): torch.manual_seed(seed) np.random.seed(seed) torch.backends.cudnn.deterministic True如果你的数据量不足十年与其硬上 GCN-LSTM不如退一步用 GCN 加轻量时序模型比如线性层加滞后特征效果往往更好也不容易出现训练不稳定的问题。6. 从复现到业务应用评估、调参与预警落地评估模型不能只看一个指标。我习惯用 RMSE 看绝对误差水平用 NSE纳什效率系数看模型相对基准的改进程度。NSE 的计算公式是 1 减去误差方差与观测方差之比NSE 大于 0.5 说明模型优于简单用历史均值预测大于 0.7 就可以用于趋势研判。业务上还要把 RMSE 和水位年变幅做对比比如年变幅 3 米的区域RMSE 在 0.3 米以内就算可用。调参不要一上来全参数网格搜索。先固定 lookback12把 KNN 的 k 扫 3、5、8 三档看验证 NSE 的差异再固定最优 k把 lookback 扫 6、12、24。小数据集上单次实验的 NSE 波动大同一组参数用三个随机种子各跑一遍取中位数再比较。这一步同时也能检验模型的稳定性——如果 k 从 3 变到 8 时 NSE 从 0.75 掉到 0.4说明模型对图结构过于敏感这个井网的边定义本身有问题需要回头检查是否有不同含水层的井被错误连边。部署层面模型输出要转成业务语言。我的做法是每月初运行一次预测输出未来 1 到 3 个月的逐井水位然后计算预测值与多年同期均值的差值即水位距平。距平低于某个阈值比如低于历史同期均值 20%就提示超采或干旱风险。把标准化统计量、邻接矩阵、井坐标和模型权重一起序列化保存部署环境换机器时直接加载不用重新计算这是我在多个项目里养成的习惯能省掉大量排错时间。这个技术方向值得投入的前提是井网有空间分布、序列长度十年以上满足这两个条件GCN-LSTM 会比单井模型带来稳定的精度提升。希望帮到你。本文还有配套的精品资源点击获取