简介这是一份面向水文地质、环境工程及给排水专业学生的毕业论文参考资料围绕《环境影响评价技术导则—地下水环境》HJ 610-2016实施中预测结果难以验证、数值法建模资料需求大的问题展开。作者以MATLAB为开发平台结合多孔介质污染物迁移动力学与MATLAB/GUI编程构建了涵盖一维、二维、三维点源瞬时注入与持续注入共六种情形的溶质运移解析模型可视化系统可输入初始浓度、孔隙度、弥散系数等参数并直观显示浓度变化并以某电厂工业废水与生活污水两类事故工况为例完成预测分析。压缩包共1个文件为pdf格式的毕业论文全文约3.02MB含中英文摘要、模型原理、界面设计与算例结果等章节结构完整适合作为建模思路、公式推导与论文写作的参照。目前已有213人学习下载可供从事地下水环境影响评价与MATLAB建模的读者借鉴。1. 地下水溶质运移预测模型为什么绕不开MATLAB设想一个很具体的场景某场地历史上堆放过含盐废渣下游 800 米处有一口供水井环保验收要回答一个问题——五年后井里的氯离子浓度会不会超过限值。这个问题不是查表能解决的它取决于地下水流速、弥散度、吸附阻滞因子和源强衰减本质上是一个对流弥散方程的定解问题。解析解只在均质、一维、定流速这类理想条件下成立真实场地的分层介质和变边界条件逼着你上数值方法。这就是“基于 MATLAB 地下水溶质运移预测模型”这个选题的现实落点把连续的偏微分方程离散成矩阵方程用时间步迭代推出浓度场的时空演化再拿实测数据把参数反演准。选 MATLAB 而不是通用编程语言理由很实在——它的矩阵语法和有限差分的天然同构sparse稀疏矩阵加反斜杠求解器几行就能解上万个节点的隐式方程组优化工具箱和神经网络工具箱又能直接把观测数据接进来标定参数画图函数顺手就能出穿透曲线和浓度等值面。所以这套东西适合三类人做水文地质或环境评价的技术人员需要一套能改参数、能出图、能写进报告的预测工具写毕业论文的学生要一个理论链条完整、代码可复现的建模流程以及想把机理模型和数据驱动方法结合的工程师。下面按理论建模、数值实现、参数反演、验证排错的顺序把每一步落到可运行的 MATLAB 代码上。2. 对流弥散方程在MATLAB里怎么推成可解的矩阵2.1 从物理过程到对流弥散控制方程地下水中的溶质同时受两件事驱动随水流整体搬运这叫对流由浓度梯度驱动向低浓度区扩散加上介质孔隙的机械弥散合起来叫水动力弥散。再叠加介质对溶质的吸附阻滞和可能的衰减一维饱和均质条件下的控制方程写作R·∂C/∂t D·∂²C/∂x² − v·∂C/∂x − λ·R·C其中 R 是阻滞因子D 是水动力弥散系数v 是孔隙流速λ 是一阶衰减常数。二维平面情形再把 y 方向的弥散项加上并把 v 写成流速矢量。这个方程是抛物线型对流占优方程难点在于对流项和弥散项的尺度差好几个数量级——对流强、弥散弱时中心差分会产生非物理振荡必须用迎风或高阶格式处理这一点在第 5 章会专门展开。从业者常见的做法是先把它无量纲化得到 Peclet 数 Pe v·Δx/D 和 Courant 数 Cr v·Δt/Δx用这两个数判断该选哪种格式、时间步能取多大。下面先看各项在 MATLAB 里怎么表达。2.2 有限差分离散从偏微分方程到三对角矩阵把求解域 [0, L] 划成 N 个节点节点间距 Δx时间用前向差分。空间二阶导用中心差分一阶导先看中心差分的形式% 一维对流弥散方程的中心差分离散系数内部节点 i 2..N-1 % C_new(i) a*C(i-1) b*C(i) c*C(i1) % 显式格式配合迎风处理对流项 dx L/(N-1); r D*dt/(R*dx^2); % 弥散项系数 pe v*dx/D; % 网格Peclet数 if pe 2 % 网格Pe过大对流项改用迎风差分避免非物理振荡 a r v*dt/(R*dx); b 1 - 2*r - v*dt/(R*dx) - lambda*dt; c r; else a r v*dt/(2*R*dx); b 1 - 2*r - lambda*dt; c r - v*dt/(2*R*dx); end这段逻辑的核心是判断网格 Peclet 数当 Pe ≤ 2 时中心差分精度更高、无数值振荡超过 2 就必须切到迎风格式用一阶精度换取稳定性。系数 a、b、c 分别对应节点 i−1、i、i1 的权重写好后整个时间推进就变成一次三对角矩阵乘法。参数 dt 是时间步长lambda 是衰减常数二者单位要和 v、D 统一到天或秒混用是新手最常踩的坑。2.3 边界条件、初始条件与稀疏矩阵组装显式格式受稳定性限制dt 通常小得可怜工程上更常用隐式或 Crank-Nicolson把未知量放到等式左边组成线性方程组 A·C^{n1} B·C^n f。这时用 MATLAB 的稀疏矩阵是关键否则 N5000 时稠密矩阵要占 200 MB 内存还解得慢。% 组装隐式格式的稀疏系数矩阵Crank-Nicolsontheta0.5 N 500; L 100; dx L/(N-1); main (1 theta*2*r)*ones(N,1); % 主对角线 off -theta*r*ones(N-1,1); % 上下对角线 A spdiags([off main off], -1:1, N, N); % 第一类边界左端定浓度右端零通量 A(1,:) 0; A(1,1) 1; b(1) C0; A(N,N-1) 0; A(N,N) 1; b(N) b(N-1); C A\b; % 反斜杠直接解稀疏线性方程组spdiags一次生成三对角稀疏矩阵A\b会自动识别稀疏结构选择合适求解器。左端 Dirichlet 边界通过置零整行并赋值实现右端 Neumann 零通量用相邻节点回代这比在矩阵里显式写差分更省事。符号含义常用单位取值来源v孔隙流速m/d达西流速除以有效孔隙度D水动力弥散系数m²/d弥散度×流速R阻滞因子无量纲1 ρb·Kd/θλ一阶衰减常数1/d现场衰减试验θ有效孔隙度无量纲室内试验或经验值表里这几个参数决定了模型的行为v 和 D 通常靠反演R 和 λ 靠试验θ 靠经验。把它们和离散系数对应起来模型才算真正落地。3. 一维溶质运移预测模型的最小可跑实现3.1 网格、时间步与Peclet、Courant数的定法动手前先做量纲分析因为它直接决定网格和时间步怎么取。实践中我一般按下面这张表拍板网格Peclet数 PeCourant数 Cr推荐格式现象 2 1中心差分显式精度高、无振荡2 ~ 10 1迎风显式有数值弥散 10任意隐式迎风或TVD弥散误差明显任意 1隐式显式必发散定网格的实用原则是让 Pe 落在 2 附近——太密浪费算力太疏直接失真。时间步按 Cr v·dt/Δx ≤ 1 反推隐式可以放宽到 Cr 5 甚至 10但 Cr 太大时会牺牲时间精度穿透曲线会显得“台阶化”。3.2 显式迎风格式与Crank-Nicolson的代码对照先用最简的显式迎风把一维模型跑通方便验证物理直觉% 一维溶质运移显式迎风格式最小实现 N 200; L 50; dx L/(N-1); v 0.5; D 0.5; R 1; lambda 0; dt 0.5*dx/v; T 100; nt round(T/dt); C zeros(N,1); C(1:5) 100; % 初始源区浓度 r D*dt/(R*dx^2); cw v*dt/(R*dx); % 迎风对流系数 for n 1:nt Cn C; for i 2:N-1 C(i) Cn(i) r*(Cn(i1)-2*Cn(i)Cn(i-1)) ... - cw*(Cn(i)-Cn(i-1)) - lambda*dt*Cn(i); end C(1) 100; % 定浓度边界 C(N) C(N-1); % 零通量边界 end内层用逐节点循环是为了看清格式结构实际推荐用矩阵向量化C Cn r*(circshift) ...速度能快十倍以上。迎风格式的特点是稳定但带数值弥散浓度锋面会被“抹平”跑出来的穿透曲线比解析解宽这是格式本身的性质不是代码错。对照实现用 Crank-Nicolson把它和显式结果叠在一张图上% Crank-Nicolson 隐式推进 theta 0.5; C zeros(N,1); C(1:5) 100; for n 1:nt % 组装右端向量依赖上一时刻C b (1-theta)*r*[0;C(1:N-2)] (1-2*(1-theta)*r)*C ... (1-theta)*r*[C(3:N);0] theta*cw*[0;C(1:N-2)] ... - (1-theta)*cw*C; b(1) 100; b(N) b(N-1); C A\b; % A为2.3节组装的稀疏矩阵 end两套代码共用同一组参数隐式可以取大十倍的时间步还稳定但单步计算量更大。判断该用哪套的标准很简单如果只是为了快速出趋势图显式够用如果要做长时间预测几年尺度隐式几乎是唯一选择。3.3 用MATLAB画图输出穿透曲线与浓度等值面模型跑完不出图等于没跑。穿透曲线看的是某个观测点浓度随时间的变化是评价报告里最核心的图件% 绘制观测点穿透曲线 x_obs round(0.6*N); % 取60%位置作为观测点 plot(t, C_hist(x_obs,:), b-, LineWidth, 1.5); hold on; plot(t, C_analytic(x_obs,:), r--); % 与解析解对照 xlabel(时间 (d)); ylabel(相对浓度 C/C0); legend(数值解,解析解); grid on; % 导出为矢量图方便插进论文 exportgraphics(gcf, breakthrough.eps, ContentType,vector);二维情形把每次时间步的浓度场用contourf或imagesc画出来再用循环生成动画或多帧图就是浓度等值面演化。配色用colormap(jet)或parula都行重点是坐标轴要标出距离和方向否则图件没有物理意义。热词里的 matlab画图、导出eps 在这一步都用得上exportgraphics比老的print更可控。4. 参数反演把BP神经网络和优化工具箱接进运移模型4.1 机理模型缺什么数据驱动补什么机理模型最大的痛点不是方程不会解而是参数不准。弥散度 α 在现场尺度上能差一到两个数量级靠经验公式给的 D 直接决定预测结果的可信度。传统做法是拿几组实测穿透曲线做最小二乘拟合但如果观测点少、地层非均质强单一套参数拟合效果会崩。数据驱动的补位思路有两种一是用 BP 网络学“观测数据到参数”的映射把反演变成回归问题二是让神经网络在机理模型基础上预测残差做混合建模。前者实现简单、易出图适合毕业论文的篇幅后者精度更高但需要更多样本。我一般先上第一种把机理模型当正演器神经网络当反演器。4.2 mapminmax归一化与BP网络拟合弥散度神经网络的输入是穿透曲线的若干采样点输出是弥散度或阻滞因子。这里必须做归一化否则输入量纲差异会让训练发散mapminmax把数据线性映射到 [−1, 1]% 用BP神经网络反演弥散度 load curves.mat % 输入: 若干穿透曲线采样点 load params.mat % 输出: 对应的弥散度标签 [in_n, in_ps] mapminmax(Curves, -1, 1); % 输入归一化 [out_n, out_ps] mapminmax(Params, -1, 1); % 输出归一化 net feedforwardnet([12 8]); % 双隐层节点数按样本量调 net.trainParam.epochs 2000; net.trainParam.lr 0.01; net.trainParam.goal 1e-4; net train(net, in_n, out_n); pred_n net(in_n); pred mapminmax(reverse, pred_n, out_ps); % 反归一化回物理量 corr(pred(:), Params(:)) % 看相关系数评估拟合质量[12 8]是隐藏层节点数样本少时取小一点否则过拟合lr学习率设在 0.01 附近较稳训练曲线一直震荡就降到 0.005。mapminmax的 ps 结构体一定要留着预测新样本时必须用同一套归一化参数否则量纲对不上结果全是错的。4.3 用lsqcurvefit和fmincon反演运移参数神经网络给的是初值真正精细标定还得回到优化工具箱。目标函数是让模型计算的穿透曲线和实测值残差平方和最小% 用lsqcurvefit反演D和R x0 [0.5, 1.2]; % 初值[D, R] lb [0.1, 0.5]; ub [5, 5]; % 参数上下界 obj (x) sim_breakthrough(x, t, dx, N, v); % 调用前向模型 [x_opt, resnorm] lsqcurvefit(obj, x0, t, C_obs, lb, ub); % 若有等式约束如质量守恒改用fmincon Aeq []; beq []; A []; b []; x_opt2 fmincon((x) sum((sim_breakthrough(x,t,dx,N,v)-C_obs).^2), ... x0, A, b, Aeq, beq, lb, ub);lsqcurvefit专门为曲线拟合设计比通用非线性求解器收敛快适合只有边界约束的场景。参数含义上x 是待反演向量lb 和 ub 卡死了物理合理范围resnorm是残差平方和用它和观测误差比较能判断拟合是否过参数化。如果还有守恒约束或参数间关系式切到fmincon把约束写进 Aeq、beq 即可。这一步和上面的神经网络形成两级反演网络给初值优化做精修初值不设好优化器很容易掉进局部极小。5. 模型精度验证与数值震荡排查5.1 质量守恒检验与解析解对照模型跑出来第一件事不是看图好不好看而是查质量守恒。数值格式只要离散方式一致总溶质质量随时间的变化应该只受衰减项影响。检验方法是逐步累加浓度场对空间积分和理论衰减曲线比% 质量守恒检验 mass_num sum(C)*dx; % 数值总质量 mass_ana mass0 * exp(-lambda*t); % 理论衰减 err abs(mass_num - mass_ana)/mass_ana; if err 0.01 warning(质量守恒偏差 %。2f%%检查边界通量或时间步, err*100); end偏差超过 1% 通常是右端零通量边界写错或者时间步太大导致显式格式的截断误差累积。另一个可靠的验证手段是拿 Ogata-Banks 一维解析解对照它描述的正是半无限域上瞬时源的浓度分布把数值解和解析解画在一张图上峰值位置对不齐说明 v 输错锋面展宽程度对不上说明 D 或格式有问题。5.2 数值震荡与数值弥散的识别和抑制两类经典误差要区分开。数值震荡表现为浓度锋面附近出现负值或超出源浓度的“波纹”根源是 Pe 2 时中心差分对流项失稳解决方式是切迎风、加密网格或换 TVD 限制器。数值弥散则相反表现为锋面被过度抹平迎风格式的一阶精度就是主因可以换二阶迎风、QUICK 格式或减小时间步缓解。排错顺序建议这样走先打印 Pe 和 Cr 确认格式选对再把网格加密一倍看结果是否收敛如果加密后解变化明显说明还没到网格无关最后检查边界处理 Dirichlet 边界 和零通量边界混用是质量不守恒的高发区。还有一个容易忽略的点是单位和量纲v 用了 m/s 而 D 用了 m²/d 是最常见的低级错误统一到天之后很多“震荡”会自己消失。模型收敛后把反演参数代回做情景预测比如源强削减 30% 后下游峰值浓度降多少这才是整套流程的最终产出。图件、参数表、守恒检验记录一起放进报告模型的可信度才算立得住。本文还有配套的精品资源点击获取
MATLAB地下水溶质运移预测模型：有限差分与参数反演
简介这是一份面向水文地质、环境工程及给排水专业学生的毕业论文参考资料围绕《环境影响评价技术导则—地下水环境》HJ 610-2016实施中预测结果难以验证、数值法建模资料需求大的问题展开。作者以MATLAB为开发平台结合多孔介质污染物迁移动力学与MATLAB/GUI编程构建了涵盖一维、二维、三维点源瞬时注入与持续注入共六种情形的溶质运移解析模型可视化系统可输入初始浓度、孔隙度、弥散系数等参数并直观显示浓度变化并以某电厂工业废水与生活污水两类事故工况为例完成预测分析。压缩包共1个文件为pdf格式的毕业论文全文约3.02MB含中英文摘要、模型原理、界面设计与算例结果等章节结构完整适合作为建模思路、公式推导与论文写作的参照。目前已有213人学习下载可供从事地下水环境影响评价与MATLAB建模的读者借鉴。1. 地下水溶质运移预测模型为什么绕不开MATLAB设想一个很具体的场景某场地历史上堆放过含盐废渣下游 800 米处有一口供水井环保验收要回答一个问题——五年后井里的氯离子浓度会不会超过限值。这个问题不是查表能解决的它取决于地下水流速、弥散度、吸附阻滞因子和源强衰减本质上是一个对流弥散方程的定解问题。解析解只在均质、一维、定流速这类理想条件下成立真实场地的分层介质和变边界条件逼着你上数值方法。这就是“基于 MATLAB 地下水溶质运移预测模型”这个选题的现实落点把连续的偏微分方程离散成矩阵方程用时间步迭代推出浓度场的时空演化再拿实测数据把参数反演准。选 MATLAB 而不是通用编程语言理由很实在——它的矩阵语法和有限差分的天然同构sparse稀疏矩阵加反斜杠求解器几行就能解上万个节点的隐式方程组优化工具箱和神经网络工具箱又能直接把观测数据接进来标定参数画图函数顺手就能出穿透曲线和浓度等值面。所以这套东西适合三类人做水文地质或环境评价的技术人员需要一套能改参数、能出图、能写进报告的预测工具写毕业论文的学生要一个理论链条完整、代码可复现的建模流程以及想把机理模型和数据驱动方法结合的工程师。下面按理论建模、数值实现、参数反演、验证排错的顺序把每一步落到可运行的 MATLAB 代码上。2. 对流弥散方程在MATLAB里怎么推成可解的矩阵2.1 从物理过程到对流弥散控制方程地下水中的溶质同时受两件事驱动随水流整体搬运这叫对流由浓度梯度驱动向低浓度区扩散加上介质孔隙的机械弥散合起来叫水动力弥散。再叠加介质对溶质的吸附阻滞和可能的衰减一维饱和均质条件下的控制方程写作R·∂C/∂t D·∂²C/∂x² − v·∂C/∂x − λ·R·C其中 R 是阻滞因子D 是水动力弥散系数v 是孔隙流速λ 是一阶衰减常数。二维平面情形再把 y 方向的弥散项加上并把 v 写成流速矢量。这个方程是抛物线型对流占优方程难点在于对流项和弥散项的尺度差好几个数量级——对流强、弥散弱时中心差分会产生非物理振荡必须用迎风或高阶格式处理这一点在第 5 章会专门展开。从业者常见的做法是先把它无量纲化得到 Peclet 数 Pe v·Δx/D 和 Courant 数 Cr v·Δt/Δx用这两个数判断该选哪种格式、时间步能取多大。下面先看各项在 MATLAB 里怎么表达。2.2 有限差分离散从偏微分方程到三对角矩阵把求解域 [0, L] 划成 N 个节点节点间距 Δx时间用前向差分。空间二阶导用中心差分一阶导先看中心差分的形式% 一维对流弥散方程的中心差分离散系数内部节点 i 2..N-1 % C_new(i) a*C(i-1) b*C(i) c*C(i1) % 显式格式配合迎风处理对流项 dx L/(N-1); r D*dt/(R*dx^2); % 弥散项系数 pe v*dx/D; % 网格Peclet数 if pe 2 % 网格Pe过大对流项改用迎风差分避免非物理振荡 a r v*dt/(R*dx); b 1 - 2*r - v*dt/(R*dx) - lambda*dt; c r; else a r v*dt/(2*R*dx); b 1 - 2*r - lambda*dt; c r - v*dt/(2*R*dx); end这段逻辑的核心是判断网格 Peclet 数当 Pe ≤ 2 时中心差分精度更高、无数值振荡超过 2 就必须切到迎风格式用一阶精度换取稳定性。系数 a、b、c 分别对应节点 i−1、i、i1 的权重写好后整个时间推进就变成一次三对角矩阵乘法。参数 dt 是时间步长lambda 是衰减常数二者单位要和 v、D 统一到天或秒混用是新手最常踩的坑。2.3 边界条件、初始条件与稀疏矩阵组装显式格式受稳定性限制dt 通常小得可怜工程上更常用隐式或 Crank-Nicolson把未知量放到等式左边组成线性方程组 A·C^{n1} B·C^n f。这时用 MATLAB 的稀疏矩阵是关键否则 N5000 时稠密矩阵要占 200 MB 内存还解得慢。% 组装隐式格式的稀疏系数矩阵Crank-Nicolsontheta0.5 N 500; L 100; dx L/(N-1); main (1 theta*2*r)*ones(N,1); % 主对角线 off -theta*r*ones(N-1,1); % 上下对角线 A spdiags([off main off], -1:1, N, N); % 第一类边界左端定浓度右端零通量 A(1,:) 0; A(1,1) 1; b(1) C0; A(N,N-1) 0; A(N,N) 1; b(N) b(N-1); C A\b; % 反斜杠直接解稀疏线性方程组spdiags一次生成三对角稀疏矩阵A\b会自动识别稀疏结构选择合适求解器。左端 Dirichlet 边界通过置零整行并赋值实现右端 Neumann 零通量用相邻节点回代这比在矩阵里显式写差分更省事。符号含义常用单位取值来源v孔隙流速m/d达西流速除以有效孔隙度D水动力弥散系数m²/d弥散度×流速R阻滞因子无量纲1 ρb·Kd/θλ一阶衰减常数1/d现场衰减试验θ有效孔隙度无量纲室内试验或经验值表里这几个参数决定了模型的行为v 和 D 通常靠反演R 和 λ 靠试验θ 靠经验。把它们和离散系数对应起来模型才算真正落地。3. 一维溶质运移预测模型的最小可跑实现3.1 网格、时间步与Peclet、Courant数的定法动手前先做量纲分析因为它直接决定网格和时间步怎么取。实践中我一般按下面这张表拍板网格Peclet数 PeCourant数 Cr推荐格式现象 2 1中心差分显式精度高、无振荡2 ~ 10 1迎风显式有数值弥散 10任意隐式迎风或TVD弥散误差明显任意 1隐式显式必发散定网格的实用原则是让 Pe 落在 2 附近——太密浪费算力太疏直接失真。时间步按 Cr v·dt/Δx ≤ 1 反推隐式可以放宽到 Cr 5 甚至 10但 Cr 太大时会牺牲时间精度穿透曲线会显得“台阶化”。3.2 显式迎风格式与Crank-Nicolson的代码对照先用最简的显式迎风把一维模型跑通方便验证物理直觉% 一维溶质运移显式迎风格式最小实现 N 200; L 50; dx L/(N-1); v 0.5; D 0.5; R 1; lambda 0; dt 0.5*dx/v; T 100; nt round(T/dt); C zeros(N,1); C(1:5) 100; % 初始源区浓度 r D*dt/(R*dx^2); cw v*dt/(R*dx); % 迎风对流系数 for n 1:nt Cn C; for i 2:N-1 C(i) Cn(i) r*(Cn(i1)-2*Cn(i)Cn(i-1)) ... - cw*(Cn(i)-Cn(i-1)) - lambda*dt*Cn(i); end C(1) 100; % 定浓度边界 C(N) C(N-1); % 零通量边界 end内层用逐节点循环是为了看清格式结构实际推荐用矩阵向量化C Cn r*(circshift) ...速度能快十倍以上。迎风格式的特点是稳定但带数值弥散浓度锋面会被“抹平”跑出来的穿透曲线比解析解宽这是格式本身的性质不是代码错。对照实现用 Crank-Nicolson把它和显式结果叠在一张图上% Crank-Nicolson 隐式推进 theta 0.5; C zeros(N,1); C(1:5) 100; for n 1:nt % 组装右端向量依赖上一时刻C b (1-theta)*r*[0;C(1:N-2)] (1-2*(1-theta)*r)*C ... (1-theta)*r*[C(3:N);0] theta*cw*[0;C(1:N-2)] ... - (1-theta)*cw*C; b(1) 100; b(N) b(N-1); C A\b; % A为2.3节组装的稀疏矩阵 end两套代码共用同一组参数隐式可以取大十倍的时间步还稳定但单步计算量更大。判断该用哪套的标准很简单如果只是为了快速出趋势图显式够用如果要做长时间预测几年尺度隐式几乎是唯一选择。3.3 用MATLAB画图输出穿透曲线与浓度等值面模型跑完不出图等于没跑。穿透曲线看的是某个观测点浓度随时间的变化是评价报告里最核心的图件% 绘制观测点穿透曲线 x_obs round(0.6*N); % 取60%位置作为观测点 plot(t, C_hist(x_obs,:), b-, LineWidth, 1.5); hold on; plot(t, C_analytic(x_obs,:), r--); % 与解析解对照 xlabel(时间 (d)); ylabel(相对浓度 C/C0); legend(数值解,解析解); grid on; % 导出为矢量图方便插进论文 exportgraphics(gcf, breakthrough.eps, ContentType,vector);二维情形把每次时间步的浓度场用contourf或imagesc画出来再用循环生成动画或多帧图就是浓度等值面演化。配色用colormap(jet)或parula都行重点是坐标轴要标出距离和方向否则图件没有物理意义。热词里的 matlab画图、导出eps 在这一步都用得上exportgraphics比老的print更可控。4. 参数反演把BP神经网络和优化工具箱接进运移模型4.1 机理模型缺什么数据驱动补什么机理模型最大的痛点不是方程不会解而是参数不准。弥散度 α 在现场尺度上能差一到两个数量级靠经验公式给的 D 直接决定预测结果的可信度。传统做法是拿几组实测穿透曲线做最小二乘拟合但如果观测点少、地层非均质强单一套参数拟合效果会崩。数据驱动的补位思路有两种一是用 BP 网络学“观测数据到参数”的映射把反演变成回归问题二是让神经网络在机理模型基础上预测残差做混合建模。前者实现简单、易出图适合毕业论文的篇幅后者精度更高但需要更多样本。我一般先上第一种把机理模型当正演器神经网络当反演器。4.2 mapminmax归一化与BP网络拟合弥散度神经网络的输入是穿透曲线的若干采样点输出是弥散度或阻滞因子。这里必须做归一化否则输入量纲差异会让训练发散mapminmax把数据线性映射到 [−1, 1]% 用BP神经网络反演弥散度 load curves.mat % 输入: 若干穿透曲线采样点 load params.mat % 输出: 对应的弥散度标签 [in_n, in_ps] mapminmax(Curves, -1, 1); % 输入归一化 [out_n, out_ps] mapminmax(Params, -1, 1); % 输出归一化 net feedforwardnet([12 8]); % 双隐层节点数按样本量调 net.trainParam.epochs 2000; net.trainParam.lr 0.01; net.trainParam.goal 1e-4; net train(net, in_n, out_n); pred_n net(in_n); pred mapminmax(reverse, pred_n, out_ps); % 反归一化回物理量 corr(pred(:), Params(:)) % 看相关系数评估拟合质量[12 8]是隐藏层节点数样本少时取小一点否则过拟合lr学习率设在 0.01 附近较稳训练曲线一直震荡就降到 0.005。mapminmax的 ps 结构体一定要留着预测新样本时必须用同一套归一化参数否则量纲对不上结果全是错的。4.3 用lsqcurvefit和fmincon反演运移参数神经网络给的是初值真正精细标定还得回到优化工具箱。目标函数是让模型计算的穿透曲线和实测值残差平方和最小% 用lsqcurvefit反演D和R x0 [0.5, 1.2]; % 初值[D, R] lb [0.1, 0.5]; ub [5, 5]; % 参数上下界 obj (x) sim_breakthrough(x, t, dx, N, v); % 调用前向模型 [x_opt, resnorm] lsqcurvefit(obj, x0, t, C_obs, lb, ub); % 若有等式约束如质量守恒改用fmincon Aeq []; beq []; A []; b []; x_opt2 fmincon((x) sum((sim_breakthrough(x,t,dx,N,v)-C_obs).^2), ... x0, A, b, Aeq, beq, lb, ub);lsqcurvefit专门为曲线拟合设计比通用非线性求解器收敛快适合只有边界约束的场景。参数含义上x 是待反演向量lb 和 ub 卡死了物理合理范围resnorm是残差平方和用它和观测误差比较能判断拟合是否过参数化。如果还有守恒约束或参数间关系式切到fmincon把约束写进 Aeq、beq 即可。这一步和上面的神经网络形成两级反演网络给初值优化做精修初值不设好优化器很容易掉进局部极小。5. 模型精度验证与数值震荡排查5.1 质量守恒检验与解析解对照模型跑出来第一件事不是看图好不好看而是查质量守恒。数值格式只要离散方式一致总溶质质量随时间的变化应该只受衰减项影响。检验方法是逐步累加浓度场对空间积分和理论衰减曲线比% 质量守恒检验 mass_num sum(C)*dx; % 数值总质量 mass_ana mass0 * exp(-lambda*t); % 理论衰减 err abs(mass_num - mass_ana)/mass_ana; if err 0.01 warning(质量守恒偏差 %。2f%%检查边界通量或时间步, err*100); end偏差超过 1% 通常是右端零通量边界写错或者时间步太大导致显式格式的截断误差累积。另一个可靠的验证手段是拿 Ogata-Banks 一维解析解对照它描述的正是半无限域上瞬时源的浓度分布把数值解和解析解画在一张图上峰值位置对不齐说明 v 输错锋面展宽程度对不上说明 D 或格式有问题。5.2 数值震荡与数值弥散的识别和抑制两类经典误差要区分开。数值震荡表现为浓度锋面附近出现负值或超出源浓度的“波纹”根源是 Pe 2 时中心差分对流项失稳解决方式是切迎风、加密网格或换 TVD 限制器。数值弥散则相反表现为锋面被过度抹平迎风格式的一阶精度就是主因可以换二阶迎风、QUICK 格式或减小时间步缓解。排错顺序建议这样走先打印 Pe 和 Cr 确认格式选对再把网格加密一倍看结果是否收敛如果加密后解变化明显说明还没到网格无关最后检查边界处理 Dirichlet 边界 和零通量边界混用是质量不守恒的高发区。还有一个容易忽略的点是单位和量纲v 用了 m/s 而 D 用了 m²/d 是最常见的低级错误统一到天之后很多“震荡”会自己消失。模型收敛后把反演参数代回做情景预测比如源强削减 30% 后下游峰值浓度降多少这才是整套流程的最终产出。图件、参数表、守恒检验记录一起放进报告模型的可信度才算立得住。本文还有配套的精品资源点击获取