☰
径向基神经网络预测地下水位:小样本场景下的稳定选择
2026/10/10 4:37:56 网站建设 项目流程

简介:基于径向基神经网络(RBFNN)的地下水位预测 MATLAB 代码,主要面向水资源管理者、环境工程人员以及机器学习初学者,用于结合降雨量、蒸发量、地表水交互、地质结构等历史数据构建地下水位预测模型,为水资源规划与环境保护提供决策支持。资源包为 zip 压缩格式,内部仅含 1 个 .m 脚本文件,整体仅 1KB,体量精炼但功能聚焦,适合快速阅读、修改与二次开发。已有 382 人浏览学习。脚本围绕数据预处理、RBFNN 网络构建、中心点确定、权值调整及预测误差评估等完整流程展开,有助于使用者理解径向基函数如何捕捉地下水位变化中的非线性特征,以及均方误差、平均绝对误差、决定系数等评估指标的应用,也可直接作为实验模板迁移至其他时间序列预测场景。借助该代码,读者能够快速上手神经网络预测建模,为地下水过度开采预警、地面沉降防治以及水资源可持续利用等实际问题提供科学参考。

1. 径向基神经网络预测地下水位:样本越少,这个老方法越值得试

手里只有一口监测井、几年逐日水位记录,领导却让你给出未来一个月的预报曲线。这种场景下,把数据直接扔给LSTM往往喂不动——样本量太少、训练不稳定。径向基神经网络(RBF)预测地下水位这条路,反而能稳定出结果。RBF网络结构简单,隐含层用一组高斯径向基函数做局部响应,中心由聚类确定、输出权重用最小二乘一步求解,训练中没有BP那种反复迭代的“玄学”。对资料有限、尺度不大的监测井项目,RBF在训练效率和泛化稳定性上常常优于BP和深网。这篇文章就按实际落地的顺序,把从数据清洗到多步预测的完整流程和参数边界讲清楚,适合水资源管理、水文地质和农业灌溉的一线数据人员照着复现。

2. RBF网络凭什么能预测地下水位:三层结构与建模选型

2.1 三层网络和“局部响应”机制

径向基神经网络的结构很直观:输入层、隐含层、输出层。隐含层的每个节点是一个径向基函数,最常用的是高斯函数 φ(||x-c||) = exp(-||x-c||² / (2σ²))。每个基函数只在中心 c 附近有显著响应,离中心越远,输出越接近 0。这种“局部响应”机制和BP网络的Sigmoid全局激励有本质区别——RBF的每个隐含节点只对输入空间的一小片区域负责,所以函数拟合的可控性更好,不容易出现BP里那种多个节点互相牵扯、训练过程震荡的情况。

输出层做的事很简单:对隐含层输出做线性加权,y = Σ w_j φ(||x - c_j||)。所以整个网络被天然拆成两部分:非线性部分由基函数的中心和宽度决定,线性部分由输出权重决定。RBF之所以训练快,就是因为它把非线性参数寻优和线性参数求解拆开了。常见的做法是第一步用K均值聚类找到中心的位置,第二步用最小二乘法直接解出权重。整个训练过程只需要一次矩阵运算,不需要像BP那样反复反向传播。

对地下水位预测这种任务,输入特征往往包含前期水位、降雨、蒸发等多个变量,特征维度不高但各变量之间存在明显的非线性耦合关系。RBF的局部响应特性正好能捕捉到“某个水位区间内,降雨响应敏感”这类局部模式。比如枯水期水位对降雨不敏感,但丰水期一场大雨水位就快速抬升,这种非线性段在不同区域表现不相同,用全局函数去拟合很容易顾此失彼,RBF靠多个局部基函数分段覆盖,反而更贴合。

2.2 中小样本场景下,为什么RBF比BP和LSTM更稳

LSTM的强项是长序列依赖,但这是用大量训练数据换来的。地下水监测井的逐日数据,实际工程里往往只有三到五年,中间还有缺测、仪器故障造成的断档。把几百条样本喂给LSTM,训练集根本撑不起门控单元的参数量,结果就是反复过拟合、验证集误差忽高忽低。BP网络在样本量小的时候也容易陷入局部极小值,十次训练十次结果不同。

RBF在这个场景下的优势是结构简单带来的稳定性。隐含节点数就是聚类中心数,一般取5到20个,参数量比LSTM小一个量级。K均值聚类给定中心之后,权重求解是凸优化问题,有唯一解,不存在局部极小。换句话讲,同一份数据用RBF训练两次,只要聚类初始化固定,结果几乎一致。这一点在实际业务里非常重要——领导让你把昨天的预测重新跑一遍,结果不能跟之前差出半米。

另一个现实因素是训练成本。地下水预测模型经常要按季度滚动重训,每次新数据进来都要重新拟合。RBF在几百条样本上训练时间是秒级和分钟级的区别,不需要GPU也能跑。这个特性让RBF特别适合嵌入到水资源管理的日常报表流程里,或者放进Excel、定时任务里定期出预报。

2.3 把水位预测建模成回归任务:输入与输出的定义

地下水位预测在数学上是一个监督回归问题。输入是当前时刻已知的变量集合,输出是未来某个时刻的地下水位埋深或标高。实际建模时,输入的构造通常遵循两个原则:物理上要和水位变化有因果关系,时序上要在预测时刻之前可获取。

最常用的输入特征包括三组:第一组是滞后水位,即前1期、前2期甚至前12期的水位观测值,用来描述系统状态惯性;第二组是气象驱动,主要是降水量、蒸发量、气温,降水是地下水补给的主要来源;第三组是人为扰动,比如农业灌溉取水量、开采井的开采量。如果所在区域有河流,还可以加上河道水位或河流流量。用公式表达就是 ŷ_{t+h} = f(x_t, x_{t-1}, ..., x_{t-n+1}),其中 n 是滑窗长度,h 是预测步长。

时间尺度的选择也很关键。逐日预测对RBF来说噪声太大,水位日变化受短时降雨影响剧烈,模型很难学到稳定规律。工程上做水资源调度,一般用旬或月尺度,滑窗长度取12个月或36旬,这样既保留了季节特征,又把高频噪声滤掉了。输出可以是单一未来值(一步预测),也可以同时输出未来多期水位(多步预测),后者在业务上更实用,但误差处理方式不同,这个放到第6章专门说。

3. 从监测井原始数据到训练样本:清洗、滑窗与特征拼接

3.1 缺测与异常水位处理:插值不能乱来

监测井数据最让人头疼的就是缺测。地下水位观测井经常因为仪器断电、传感器故障、人为疏忽出现连续几天甚至几周的空缺。处理缺测的第一步是先看缺测的分布——如果缺测集中在某个季节,线性插值会把那个季节的真实波动抹平,所以插值前要先按月份统计缺失比例。

我一般用的是两层处理:短缺测(7天以内)用线性插值,因为地下水位短期变化是连续过程,线性插值带来的误差在可接受范围;长缺测(超过30天)直接放弃该段,不参与样本构造,避免插值引入的虚假趋势。异常值的识别用相邻水位差法——地下水位自然变化在非强降雨期每天不会超过0.3到0.5米,如果相邻两日水位差超过这个阈值,优先怀疑是仪器跳变或人为录入错误。

import pandas as pd import numpy as np # 读取监测井逐日水位数据 df = pd.read_csv('well_001.csv', parse_dates=['date'], index_col='date') df = df[['水位']].copy() # 短缺测线性插值:limit=7表示最多连续插7天,超过则保留缺失 df['水位'] = df['水位'].interpolate(method='linear', limit=7) # 异常值检测:相邻水位差超过0.5米视为仪器跳变 diff = df['水位'].diff().abs() bad = diff > 0.5 df.loc[bad, '水位'] = np.nan # 二次插值回填,过滤掉跳变点 df['水位'] = df['水位'].interpolate(method='linear', limit=7)

注意这里两个参数需要根据研究区实际情况调整。0.5米阈值适用于平原区潜水井,如果在山区基岩裂隙水区域,水位对降雨响应快,日变化可能超过1米,阈值要放大。limit=7的选择逻辑是:地下水位的自相关系数在7天尺度上通常还在0.8以上,线性插值可信度较高;超过7天,宁可让样本缺失,也不要制造一段假数据。

插值完成后还有一个容易忽略的步骤:检查清洗后的水位序列和降雨记录的对应关系。一个简单的验证方法是,把水位上升段和前期降雨日对齐,如果存在“连续大雨后水位完全不动”的片段,那段时间的数据大概率有问题,需要回到原始台账核对。

3.2 滑动窗口构造样本:时间步怎么定、留多少预测步

清洗完原始序列后,下一步是把时间序列转成监督学习的样本对。核心是滑窗函数:用过去 n 个时刻的特征预测未来 h 时刻的水位。窗口长度 n 的选择直接决定模型能看到多大的历史信息。月尺度数据我一般取 n=12,对应12个月的季节性周期;旬尺度数据取 n=36,同样覆盖一整年。

构造样本时必须注意两个细节。第一,滑窗步长通常是1个时间单位,也就是逐月滑动,这样样本量等于“总时序长度 - n - h + 1”,几百条水位记录能产生几百个样本;第二,样本之间是有重叠的,这会导致训练集内部存在信息冗余,但不影响模型训练,只是评估时要按时间块切分,这个在第4章专门说。

def make_sequences(data, feature_cols, n_in=12, n_out=1): """ 将时序数据转为滑窗样本 data: DataFrame,包含水位和气象特征 feature_cols: 参与建模的特征列名 n_in: 输入窗口长度(月/旬) n_out: 输出步长(预测未来n_out期) """ X, y = [], [] for i in range(len(data) - n_in - n_out + 1): # 取i到i+n_in时刻的特征作为输入 X.append(data.iloc[i:i+n_in][feature_cols].values) # 取i+n_in到i+n_in+n_out时刻的水位作为输出 y.append(data.iloc[i+n_in:i+n_in+n_out]['水位'].values) return np.array(X), np.array(y) # 特征列:水位、降雨、蒸发 feature_cols = ['水位', '降雨量', '蒸发量'] X, y = make_sequences(df, feature_cols, n_in=12, n_out=1) print(X.shape, y.shape) # (样本数, 12, 3) (样本数, 1)

这里有几个工程参数需要注意。n_in=12是默认起点,但如果你的监测井位于强开采区,水位受开采影响远大于气候影响,那么半年内的开采量数据比一年前的水位更重要,此时n_in=6可能更好。n_out的取值决定了模型的预测目标:n_out=1做单步预测,误差最小;n_out=3或6做月尺度的多步预测,业务上更实用,但误差会随步长增加。特征列的拼接顺序保持一致,预测时你提供什么特征顺序,模型就按什么顺序读取。

3.3 时序分解给RBF“减负”:先拆掉季节项,再学残差

地下水位序列最大的特点是强季节性——北方平原区春季农业开采导致水位下降,夏季雨季补给后水位回升,这种周期波动幅度经常超过总方差的60%。如果让RBF直接拟合原始水位,基函数的大部分能力都花在学习这条固定的季节性曲线上了,剩下对降雨、开采等异常事件的响应能力就变弱了。

一个有效的做法是参考时序分解的思想,把水位序列拆成趋势项、季节项和残差项,让RBF只预测残差部分。这样基函数不需要用大量节点去逼近一条平滑的季节曲线,可以把表达力集中在“偏离常态”的部分——而这恰恰是预测真正有价值的地方。在实际项目中,我用STL分解后,预测误差能下降20%到30%。

分解后的预测策略是:季节项和趋势项外推(季节项按周期延拓、趋势项用线性回归外推),残差项用RBF预测,最后把三部分加起来得到最终水位。这个流程下的RBF实际上成了一个残差回归器,任务从“预测一条大波动曲线”简化为“预测偏离量”,拟合难度显著降低。如果不想引入STL依赖,也可以使用季节性差分(即当前水位减去一年前同期水位),这样也能剥离掉大部分季节性,操作更简单。

4. 用Python实现RBF水位预测模型:聚类定中心、伪逆定权重

4.1 最小可用的RBF网络实现:不到60行代码

现在到核心环节——手写一个可用的RBF网络。我这里用的是numpy手动实现,不放scikit-learn里的SVR带RBF核去替代,原因是SVR的求解目标是最小化ε不敏感损失,和RBF网络的直接回归在数学上不是一回事,工程上手动实现反而可控性更强。

实现分三步:第一步用K均值聚类确定基函数中心;第二步根据中心之间的距离计算基函数宽度σ,这里用中心距离中位数乘以spread系数来得到;第三步构造隐含层输出矩阵H,用伪逆求解输出权重。伪逆求解对比直接求逆的好处是,当基函数输出矩阵病态时(比如两个中心距离太近导致高度相关),伪逆仍然能给出有界的最小二乘解。

import numpy as np from sklearn.cluster import KMeans class RBFNetwork: def __init__(self, n_centers=10, spread=0.8): self.n_centers = n_centers self.spread = spread self.centers = None self.weights = None self.sigma = None def _basis(self, X): # 计算高斯径向基输出矩阵 # X: (n_samples, n_features) -> H: (n_samples, n_centers) diff = X[:, None, :] - self.centers[None, :, :] d2 = np.sum(diff ** 2, axis=-1) return np.exp(-d2 / (2 * self.sigma ** 2)) def fit(self, X, y): # 第一步:K均值聚类找中心 km = KMeans(n_clusters=self.n_centers, random_state=0, n_init=10) km.fit(X) self.centers = km.cluster_centers_ # 第二步:用中心间距离的中位数定sigma,spread是缩放系数 dists = np.linalg.norm(self.centers[:, None, :] - self.centers[None, :, :], axis=-1) median_dist = np.median(dists[dists > 0]) self.sigma = self.spread * median_dist # 第三步:构造隐层输出矩阵,伪逆求解输出权重 H = self._basis(X) self.weights = np.linalg.pinv(H) @ y def predict(self, X): H = self._basis(X) return H @ self.weights

这段代码有几个关键设计。宽度σ不是用所有中心距离的平均值,而是用中位数,因为聚类后个别中心可能落在数据稀疏区,离其他中心很远,会把平均值拉大导致基函数整体变宽、局部分辨率下降。spread参数是唯一的自由超参数,它控制基函数的覆盖范围——spread越大,每个基函数的影响范围越广,函数曲线越平滑;spread越小,基函数越尖锐,拟合能力越强但越容易过拟合。n_centers控制模型容量,相当于BP网络里隐含层节点数。

还有一个容易被忽略的工程细节:滑窗构造出的X是三维张量(样本数、时间步、特征数),而RBF网络的输入层需要展平成一维向量。所以在调用fit之前必须把三维张量reshape成二维(样本数, 时间步×特征数)。这一步做漏了,代码会在距离计算时报维度错误。

# 滑窗得到的X是(样本数, 12, 3),要展平成(样本数, 36) n_samples = X.shape[0] X_flat = X.reshape(n_samples, -1)

4.2 训练与验证切分:时间序列不能用随机打乱

这是RBF水位预测项目里最常见的错误:用train_test_split随机切分训练集和验证集,结果验证集里混着训练集前后的时间点,因为滑窗重叠,模型在验证集上的表现虚高。地下水位是强自相关序列,今天的样本和30天前的样本高度相关,随机切分等于让模型提前“看到”了答案。

正确做法是按时间顺序切分:前80%的样本训练,后20%的样本验证。如果要更严谨,用TimeSeriesSplit做滚动交叉验证——每次用前面的块做训练、后面的块做验证,这样模拟的就是“用历史预测未来”的真实场景。

from sklearn.model_selection import TimeSeriesSplit def rmse(y_true, y_pred): return np.sqrt(np.mean((y_true - y_pred) ** 2)) def nse(y_true, y_pred): # Nash-Sutcliffe效率系数,1表示完美拟合,<0表示比用均值预测还差 return 1 - np.sum((y_true - y_pred) ** 2) / np.sum((y_true - np.mean(y_true)) ** 2) # 按时间顺序切分,前80%训练,后20%验证 split = int(n_samples * 0.8) X_train, X_test = X_flat[:split], X_flat[split:] y_train, y_test = y[:split], y[split:] model = RBFNetwork(n_centers=10, spread=0.8) model.fit(X_train, y_train) y_pred = model.predict(X_test) print(f"RMSE: {rmse(y_test, y_pred):.3f} m") print(f"NSE: {nse(y_test, y_pred):.3f}")

评价指标上,RMSE的单位是米,直接反映预测误差的物理量级,和领导汇报时用这个最直观。NSE是无量纲指标,含义是模型预测比“直接用历史均值预测”好多少——NSE大于0.6说明模型有使用价值,大于0.8说明精度良好。如果NSE低于0.5,先不要急着调模型参数,回头检查数据质量和特征构造。另外要打印预测值和实测值的散点对比,重点关注有没有系统性低估高峰水位的问题。

4.3 SPREAD与聚类数的调参边界:一个网格搜索搞定

RBF网络仅有的两个超参数是n_centers和spread,但这两个参数的组合对结果影响很大。工程上不需要复杂的贝叶斯优化,用网格搜索加TimeSeriesSplit交叉验证就足够。

import itertools # 参数网格 param_grid = { 'n_centers': [5, 8, 10, 15, 20], 'spread': [0.3, 0.5, 0.8, 1.2, 2.0] } best_score = float('inf') best_params = None tscv = TimeSeriesSplit(n_splits=5) for n_centers, spread in itertools.product(param_grid['n_centers'], param_grid['spread']): scores = [] for train_idx, val_idx in tscv.split(X_flat): model = RBFNetwork(n_centers=n_centers, spread=spread) model.fit(X_flat[train_idx], y_train[train_idx]) pred = model.predict(X_flat[val_idx]) scores.append(rmse(y_true[val_idx], pred)) mean_score = np.mean(scores) if mean_score < best_score: best_score = mean_score best_params = (n_centers, spread) print(f"最优参数: n_centers={best_params[0]}, spread={best_params[1]}, RMSE={best_score:.3f} m")

调参时有两条边界经验。n_centers的上限不要超过训练样本量的10%,否则每个基函数只覆盖少数几个样本,拟合的是噪声而不是规律。spread的常见有效区间是0.3到2.0,如果最优spread出现在边界值(小于0.1或大于5),说明数据分布有问题,比如特征没做归一化,或者聚类中心分布极不均匀。这里强烈建议对输入特征做MinMax归一化到0到1区间后再训练,否则降雨量(几百毫米)和水位(几十米)的量纲差异会让距离计算完全由大数值特征主导。

5. RBF水位预测避坑指南:五个真实踩坑记录

5.1 归一化泄漏:用全序列算MinMax,你已经偷看了未来

现象:验证集上RMSE只有0.05米,模型精度漂亮得惊人,但一投入实际预测就完全不准,误差超过0.5米。

原因:做MinMax归一化时,用了整个数据集(包括测试集)的最大值最小值。这等于归一化阶段就把测试集的信息泄露给了训练过程。地下水位序列有趋势性,测试集往往位于序列末端,其水位最大值普遍高于训练集,归一化后测试集的数值范围被压缩,模型预测结果自然偏向训练集均值。

解决:归一化参数只能用训练集的统计量,然后分别应用于训练集和测试集。更稳妥的做法是先把数据按时间切分,再在训练集上fit归一化器,在测试集上只做transform。如果做滚动预测,每次重训都要重新fit归一化器,不能沿用上一轮的参数。

5.2 聚类数K的“玄学”:太少欠拟合,太多过拟合

现象:n_centers从5调到15时,训练集误差持续下降,但验证集误差在10以后开始反弹,NSE从0.72掉到0.58。

原因:地下水位变化不是均匀的——平水期曲线平滑,汛期急涨急落。聚类数太少时,基函数中心分布在数据密集区,汛期水位变化大的区域没有中心覆盖,拟合不出尖峰;聚类数太多时,平水期的样本被分成很多小簇,每个基函数只匹配少量样本,学到的全是细节噪声。

解决:不要光看训练误差曲线,要同时画训练集和验证集的误差曲线,找两者的交叉点。另外观察聚类中心在特征空间里的分布——如果大量中心挤在同一个区域、个别区域完全没有中心,说明K值不合适,或者某些特征对距离计算贡献过大,需要检查特征归一化。

5.3 SPREAD越大误差越低?那是假象

现象:把spread从0.5调到2.0,训练集RMSE从0.15降到0.08,看起来拟合越来越好,但验证集RMSE从0.20涨到0.40,模型完全翻车。

原因:spread增大导致每个基函数覆盖范围变宽、函数曲线变平滑,在训练集上确实降低了残差。但spread过大时,所有基函数几乎重叠成同一个“山丘”,网络表达能力退化为一个宽平滑函数。这时候RBF没有了局部响应特性,本质上回到了全局拟合,对输入特征的细微变化完全不敏感。

解决:把spread当成正则化参数来看待,而不是拟合精度旋钮。正确的调参目标是让验证集误差最小,而不是训练集误差最小。判断spread是否过大,可以打印隐含层输出矩阵H的秩——如果H的秩远小于n_centers,说明基函数严重冗余,spread过大了。

5.4 递归多步预测误差累积:自己预测的结果喂给自己,越走越偏

现象:单步预测RMSE只有0.12米,但滚动预测未来6个月时,到第4个月水位预测曲线就成了水平直线,完全失去波动。

原因:递归多步策略是把第t+1步的预测值当作输入,去预测第t+2步。RBF本质上是回归模型,预测值会向均值回归,把预测值再喂回去,相当于输入信号本身就在被“钝化”。每走一步,信号的高频成分被削弱一层,误差逐层累积,最终输出趋于平滑直线。

解决:多步预测改用直接策略——把输出维度从1改成h,让模型一次输出未来h期的水位。代码里把make_sequences的n_out设为6,同时把RBF输出层的weights矩阵维度变成(n_centers, 6),训练时y变成(n_samples, 6)。代价是输出维度增加,验证误差会比单步高一些,但不会出现误差递归累积导致的曲线“死掉”问题。

5.5 单井模型直接搬到邻井,结果惨不忍睹

现象:A井模型训练效果良好,NSE达到0.85;直接拿到距A井仅2公里的B井上去预测,NSE直接掉到0.1以下。

原因:地下水位受局部水文地质条件控制——A井可能位于砂层富水区,降雨入渗补给快;B井可能在粘土层覆盖下,降水补给比例极低且滞后严重。两个井的水位变幅和相位差完全不是一套规律,特征和输出之间的映射关系不同。所有回归模型都只能内插,不能外推到训练数据覆盖的范围之外。

解决:要么每个井单独建模(这是最可靠的做法),要么做多井联合建模时把井的空间属性(含水层类型、距河流距离、开采强度、井深)作为特征输入进模型。空间属性是静态特征,在滑窗内逐时间步重复拼接即可。如果区域内井数量多、数据都足够,联合建模才有意义;否则老老实实一井一模型。

6. 从单步到多步:把RBF预测接到水资源业务闭环里

模型训练好、参数调优之后,真正的价值在于产出对未来水位过程线的预测,而不是停留在论文里的那个RMSE数字上。在这个环节,我常用的做法是“直接多步预测 + 月度滚动重训”。每次预测时向模型输入最近12个月的水位、降雨、蒸发数据,一次输出未来3个月的水位过程,然后每月月初用新观测数据重新训练一次。这样既避免递归预测的误差累积,又保证模型参数不落后于当前水文条件。

验证方法上,滚动原点预测比一次性切分更贴近业务:从第24个月开始,每个月生成一组未来3个月预测,把多组预测拼成一条连续曲线,再与实测值对比。计算整条曲线的NSE和峰值水位误差,峰值误差尤其重要——防洪和生态补水决策都依赖汛期高水位是否报得准。如果峰值系统性偏低,可以在模型输出后加一个经验修正项,比如研究区内降雨大于50毫米后水位上升幅度的统计值。

另一个值得养成的习惯是记录每次预测的输入数据范围和预测结果,建立“模型预测档案”。这个档案可以在模型预测出现偏差时回溯原因:是降雨特征没带上,还是监测井水位计该校准了。早期我吃过这方面的亏——模型预测连续偏低一周,后来排查发现是监测井的传感器漂移了,不是模型坏了。数据质量的问题不能靠调参数解决。地下水位预测这件事,模型只是最后一步,前面的数据功夫才是决定成败的地方。希望这些经验能帮你在自己的项目里少走几步弯路,缩短从数据到决策之间的距离。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询