简介一份面向水文专业学生、水利工程技术人员及MATLAB初学者的PDF学习资料聚焦如何用MATLAB处理水文计算中的常见问题。内容围绕单位线推求、相关分析、系列插补延长三大应用展开重点讲解最小二乘法推求单位线并将流量过程表达式转化为矩阵形式求解附带具体实例与计算表便于读者按步骤理解编程思路和数学模型转换方法。文中演示了矩阵运算在单位线求解中的具体实现对理解最小二乘估计与水文统计相结合很有帮助。全文为期刊论文格式的参考文献可作为课程作业、毕业设计或实际工程计算的参考资料。资源包共1个PDF文件大小仅153KB轻量易下载适合移动端随时阅读已有270人学习使用对于希望借助MATLAB提升水文计算效率、降低手工计算工作量的读者具有实用价值。1. 用最小二乘法反推单位线MATLAB 把水文计算从手算中解放出来水文计算里最磨人的不是公式本身而是反复解矩阵、做相关分析、插补延长系列这些体力活。拿单位线推求来说传统方法需要用手算分解地面径流过程线每来一场洪水就要重新解一遍矛盾方程组一个流域多条单位线综合下来工作量能占掉设计周期的大半。这篇 2005 年发表在水文计算领域的经典应用文章讲的是用 MATLAB 把单位线推求、相关分析、系列插补延长三件事变成十几行命令核心思路是把卷积关系写成矩阵方程 h·q Q再用最小二乘一次解出单位线纵坐标不需要写循环、不需要迭代调参数据准备好以后计算过程可以压缩到几秒。文章适合两类人一是做工程水文设计、需要频繁推求单位线和延长径流系列的水利工程师二是准备用 MATLAB 做科学计算、想找真实案例而不是玩具 demo 的初学者。这里的价值不只是省时间而是把统计规律 成因规律这两条水文分析主线用矩阵运算统一起来思路本身比命令更值得拆解。2. 单位线推求的矩阵原理与最小二乘实现2.1 从径流过程到线性方程组为什么单位线推求本质上是解卷积单位线的定义很简洁给定流域上单位时段内均匀分布的单位净雨深通常取 10mm在出口断面形成的地面径流过程线。它把净雨 → 径流的转换关系线性化了所以一旦有了单位线任意净雨过程都能通过卷积叠加推出流量过程反过来有了实测净雨和地面径流过程也能反推单位线。反推的数学形式是线性方程组。设地面流量过程有 l 个时段净雨有 m 个时段单位线有 n 个时段三者满足 n l - m 1。把卷积展开Q1 h1·q1 Q2 h2·q1 h1·q2 Q3 h3·q1 h2·q2 h1·q3 ... Ql hm·q(n)这里的 h 是净雨量q 是单位线纵坐标Q 是地面径流量。用矩阵写就是h·q Q其中 h 是 l 行 n 列的带状矩阵每一行对应一个时段的地面流量方程。工程上实测数据通常 m 大于 1方程组个数多于未知数 q 的个数也就是矛盾方程组没有精确解只能在最小二乘意义下找最优解。这就是整个推求过程的核心难点不是不会列方程而是方程列出来解不干净。2.2 正规方程与左除避免显式求逆的两种解法传统教科书给的标准做法是把矛盾方程组转成正规方程也就是用 h 的转置左乘两端hᵀ·h·q hᵀ·Q q (hᵀ·h)⁻¹·hᵀ·Q这在理论上没有问题但实际用 MATLAB 计算时我一般不会直接写 inv(h*h)*h*Q原因有两个一是 hᵀh 的条件数是 h 条件数的平方如果 h 本身病态求逆会把误差放大二是显式求逆在矩阵规模变大时浪费算力。更稳的写法是用左除q h \ Q;左除运算符 \ 会根据矩阵结构自动选择 LU 分解、Cholesky 分解或 QR 分解数值稳定性比显式求逆好得多这是 MATLAB 里最值得养成习惯的一个细节。如果确实需要看到正规方程的中间形态可以这样写A h * h; % 正规方程的系数矩阵 b h * Q; % 正规方程的右端项 q_normal A \ b; % 用左除解正规方程这里的 q_normal 和直接用 h\Q 得到的结果理论上是一致的但后者不放大条件数。实际计算时如果发现两种结果有明显差异说明 h 矩阵病态需要检查净雨时段划分是否合理、流量过程是否已经正确分割基流。2.3 实例复现构造 h 矩阵与 Q 向量论文给了一个完整实例某流域实测 20 个时段的地面径流量净雨过程有两个时段15.7mm 和 5.9mm时段长 Δt 12h要求推求单位线。按 n l - m 1单位线时段数为 20 - 2 1 19 个。构造矩阵的关键是理解 h 的结构第 1 列是净雨序列从第 1 时段开始依次下移第 2 列从第 2 时段开始下移以此类推。完整的 MATLAB 代码如下% 构造净雨矩阵 h20行 × 19列 % 第1列15.7 从第1行开始之后依次递推 % 第2列5.9 从第1行开始15.7 从第2行开始 h [15.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 5.9, 15.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 5.9, 15.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 5.9, 15.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 5.9, 15.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 5.9, 15.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 5.9, 15.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 5.9, 15.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 5.9, 15.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 5.9, 15.7, 0, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.9, 15.7, 0, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.9, 15.7, 0, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.9, 15.7, 0, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.9, 15.7, 0, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.9, 15.7, 0, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.9, 15.7, 0, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.9, 15.7, 0, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.9, 15.7, 0; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.9, 15.7; 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5.9]; % 实测地面径流过程 Qm³/s共21个值 Q [0; 20; 275; 737; 1065; 840; 575; 389; 261; 180; 128; 95; 73; 55; 40; 29; 19; 12; 6; 1; 0]; % 最小二乘推求单位线纵坐标 q h \ Q; % 单位净雨深为10mm换算为10mm的单位线 q_10mm q * 10; % 绘制单位线 figure; plot(q_10mm, o-, LineWidth, 1.5); xlabel(时段数 (Δt 12h)); ylabel(单位线纵坐标 q (m³/s)); title(最小二乘法推求单位线); grid on;代码里的 h 矩阵是 20 行 19 列第 1 列由 15.7 开头第 2 列由 5.9 开头两列之间错位一个时段这正好对应两时段净雨对出口流量的叠加关系。Q 向量是地表径流过程已经完成了地下径流分割这一步不能省否则推出来的单位线会带上基流的尾巴。q 向量是单位净雨深取 1mm 时的结果要得到标准单位线需要按实际单位净雨深缩放也就是乘 10。矩阵/向量维度含义hl × n净雨矩阵带状结构每列对应一个净雨时段的贡献Ql × 1实测地面径流过程需先分割基流qn × 1待求单位线纵坐标1mm 净雨深q_10mmn × 1换算为标准 10mm 单位线2.4 负值修正与光滑化最小二乘解不直接可用怎么办直接解出来的 q 向量很可能出现负值这在水文物理上不成立单位线纵坐标不可能为负。负值的原因是方程组矛盾程度高、净雨资料存在误差或者基流分割不够彻底。常见做法是把负值置零然后对非零段做一次光滑处理再重新用实测流量过程检验推流结果。% 负值修正将负值置零 q_corrected max(q, 0); % 三点滑动平均光滑 q_smooth conv(q_corrected, ones(1,3)/3, same); q_smooth max(q_smooth, 0); % 检验用修正后的单位线重新推流对比实测值 Q_sim conv(h * 10, q_corrected / 10); Q_sim Q_sim(1:length(Q)); % 计算纳什效率系数 NSE NSE 1 - sum((Q - Q_sim).^2) / sum((Q - mean(Q)).^2); fprintf(NSE %.3f\n, NSE);需要说明的是这里的单位换算容易出错。h 矩阵里是毫米单位的净雨Q 是立方米每秒的流量两者之间还有一个流域面积换算系数实际应用中需要乘以流域面积并做单位统一论文里为书写方便把单位净雨深取为 1mm工程计算时要还原回去。3. 相关分析与回归方程polyfit 和 corrcoef 的工程用法3.1 水文相关分析的问题背景水文现象之间的关系不是严格的函数关系比如年降雨量和年径流量同一条流域、同一量级的降雨径流量也会有波动因为还受蒸发、土壤前期含水量、下垫面条件影响。但两者确实存在统计相关性可以用回归方程描述。相关分析在水文计算里有两个典型用途一是插补缺测年份的径流资料利用长系列降雨资料延长短系列径流资料二是检验两个水文变量之间的关联强度比如上下游站的流量相关。论文里的例子就是典型的降雨-径流相关延长表 2 给了 1954 到 1965 年的年降雨量和年径流量数据要用这 12 对数据建立回归方程然后利用已知的年降雨量推算缺测年份的年径流量。3.2 直线回归的最小二乘原理与 polyfit 对应关系两变量直线相关的基本模型是 y a bx其中 x 是自变量y 是倚变量。最小二乘原理是让观测点与拟合直线在纵轴方向的离差平方和最小也就是S Σ(yi - a - b·xi)² → min对 a 和 b 分别求偏导并令其为零得到两个正规方程。解出来之后回归系数 b 可以用相关系数 r 和两个系列的均方差表示截距 a 用均值关系确定。推导结果在论文里有完整展示这里直接说结论回归方程本质上是在用 x 的变异信息去估计 y 的期望值相关系数 r 衡量的是线性关系的紧密程度。在 MATLAB 里polyfit 函数做多项式拟合n1 时就是直线拟合返回的两个值分别是斜率 b 和截距 a顺序和手动推导的公式正好反过来初学者经常在这里搞混。corrcoef 函数返回的是相关系数矩阵对角线是 1非对角线是两个变量之间的相关系数。3.3 完整代码回归、相关系数、插补延长一条龙% 年降雨量序列 x (mm) x [2014, 1211, 1728, 1157, 1257, 1029, 1306, 1029, 1310, 1356, 1266, 1052]; % 年径流量序列 y (mm) y [1362, 728, 1369, 695, 720, 534, 778, 337, 809, 929, 796, 383]; % 直线拟合polyfit(x, y, 1) 返回 [斜率 b, 截距 a] p polyfit(x, y, 1); b p(1); % 斜率 a p(2); % 截距 fprintf(回归方程: y %.4f·x %.4f\n, b, a); % 计算相关系数 R corrcoef(x, y); r R(1, 2); fprintf(相关系数 r %.4f\n, r); % 绘制散点图与拟合线 figure; scatter(x, y, 40, filled); hold on; x_fit linspace(min(x), max(x), 100); y_fit polyval(p, x_fit); plot(x_fit, y_fit, r-, LineWidth, 1.5); xlabel(年降雨量 x (mm)); ylabel(年径流量 y (mm)); title(降雨-径流相关分析); grid on; % 插补延长已知缺测年份的降雨量 xi [998; 1023; 1459; 1327.5]; yi a b * xi; % 输出插补结果 fprintf(插补延长的年径流量:\n); disp([xi, yi]);参数说明polyfit 的第三个参数是多项式阶数1 表示直线2 表示抛物线水文相关分析绝大多数场景用 1 阶就够盲目用高阶拟合容易过拟合。corrcoef(x,y) 返回 2×2 矩阵取 R(1,2) 是因为这里只关心 x 和 y 的相关系数如果传一个多列矩阵进去返回的就是两两之间的相关系数矩阵适合分析多变量相关性。论文实例的结果是回归方程 y 1.0477x - 585.3602相关系数 r 0.9503。r 的平方是 0.90 左右说明年径流量的变异大约有 90% 能被年降雨量解释这个量级的相关性在水文上算相当理想。3.4 相关系数显著性检验r 大不等于真的相关r 0.95 看起来很漂亮但如果样本太少这个相关可能是偶然的。水文上一般要求做显著性检验最常用的是 t 检验n length(x); df n - 2; % 自由度 t_stat r * sqrt(df) / sqrt(1 - r^2); p_value 2 * (1 - tcdf(abs(t_stat), df)); fprintf(t统计量 %.3f, p值 %.4f\n, t_stat, p_value); % 查表临界值显著性水平 α 0.05 t_crit tinv(0.975, df); fprintf(95%%置信水平临界值 t(0.025, %d) %.3f\n, df, t_crit);对 12 个样本、自由度 1095% 置信水平的 t 临界值大约是 2.228这里算出来的 t 统计量落在拒绝域里可以认为相关显著。p 值的意义是如果两个变量真的不相关出现当前样本相关系数的概率只有 p_value 这么多p 越小越有理由相信相关是真实的。4. 系列插补延长与结果验证从回归方程到可用设计成果4.1 插补延长的可靠边界什么情况下能外推什么情况下不能有了回归方程把缺测年份的降雨量代进去就能得到径流量。这个操作技术上很简单工程上却要谨慎。回归方程的适用区间是参与拟合的 x 数据范围也就是 1029mm 到 2014mm 之间在这个范围内内插是可靠的超出范围外推要小心因为降雨-径流关系在高雨量段可能偏离线性比如超大暴雨年份径流系数会变化。论文里的插补值包括 998mm、1023mm、1459mm 和 1327.5mm其中 998mm 和 1023mm 略低于实测样本的最小值 1029mm属于轻微外推在工程上可以接受但如果要外推到 500mm 或者 3000mm就应该谨慎或者改用非线性拟合。这是使用回归方程做插补延长时最容易被忽略的问题。4.2 相对误差与残差分析检验延长结果的合理性插补出来的值是点估计没有给出不确定性范围。工程上可以用回归的标准误差来衡量估计精度% 残差计算 y_fitted a b * x; residual y - y_fitted; % 回归标准误差 Sy Sy sqrt(sum(residual.^2) / (length(x) - 2)); % 各插补点的95%预测区间简化方式 for i 1:length(xi) se_pred Sy * sqrt(1 1/length(x) (xi(i) - mean(x))^2 / sum((x - mean(x)).^2)); half_width tinv(0.975, length(x)-2) * se_pred; fprintf(xi%.1fmm, 预测 y%.1fmm, 95%%区间 [%.1f, %.1f]\n, ... xi(i), a b*xi(i), a b*xi(i) - half_width, a b*xi(i) half_width); end这里的预测区间比单纯的回归线置信区间宽因为它包含了新观测点的随机误差。水文上常用这个概念来评估插补值的可靠程度如果延长出来的径流量落在合理范围内可以用如果区间宽到失去工程意义就要考虑增加参证站或者用其他方法交叉验证。4.3 交叉验证用留一法检验回归方程的稳定性样本量不大时我更推荐做一次留一交叉验证Leave-One-Out Cross-Validation, LOOCV每次拿掉一个样本用剩下的 11 个样本拟合然后预测被拿掉的那个点最后统计误差x x(:); y y(:); n length(x); pred zeros(n, 1); for i 1:n train_idx setdiff(1:n, i); p_train polyfit(x(train_idx), y(train_idx), 1); pred(i) polyval(p_train, x(i)); end % 计算均方根误差和平均相对误差 rmse sqrt(mean((y - pred).^2)); mape mean(abs((y - pred) ./ y)) * 100; fprintf(LOOCV RMSE %.2f mm\n, rmse); fprintf(LOOCV MAPE %.1f%%\n, mape);LOOCV 的结果能反映出回归方程对单个样本的敏感程度。如果删除某个样本后拟合结果变化很大说明这个样本是强影响点需要检查它是否有特殊的物理成因比如某一年发生了特大暴雨或干旱。在水文分析里这种稳健性检验比单看 r 值更有工程意义因为延长出来的系列会直接参与频率计算一个不稳定回归方程会给设计洪水结果带来隐性偏差。论文里没有做这一步但以现在的水文计算规范要求频率分析之前对插补延长系列做合理性检查是标准流程。把 LOOCV 作为验证步骤加进去等于给旧方法补齐了现代工程实践的可靠性要求。4.4 从 MATLAB 脚本到完整分析流程的整理把这一套在 MATLAB 环境里整合起来一个完整的水文相关分析脚本可以按这样的结构组织模块内容输出数据准备录入实测序列分割基流x、y 列向量单位线推求构造 h 矩阵最小二乘解 q单位线纵坐标相关性分析polyfit、corrcoef、显著性检验回归方程、r、p 值插补延长代入缺测年份降雨量计算预测区间延长后的径流系列验证LOOCV、残差分析、过程线对比RMSE、MAPE、合理性结论脚本里把每个模块的输出都用 fprintf 打印出来方便直接粘贴到计算书里。这也是这类计算容易被忽略的一点水文计算要求可追溯性脚本里每一步都应该有文档化输出而不是只保留最终结果。实际项目中我习惯在脚本开头加一段注释块记录资料系列、流域面积、基流分割方法这样三个月后回头检查时还能复现当时的计算条件。用 MATLAB 做这类工作真正的生产力提升不只在解方程那一下而是把整个分析流程变成一个可重复执行的脚本这也是这篇 2005 年的方法放在今天依然有参考价值的原因。本文还有配套的精品资源点击获取