简介一份完整的数学建模竞赛A题参赛文档以城市表层土壤重金属污染分析为项目背景系统讲解基于Matlab的环境数据建模与可视化流程。内容涵盖题目重述、问题分析、模型假设、模型建立与求解、模型评价与推广等完整环节重点展示As、Cd、Cr、Cu、Hg、Ni、Pb、Zn这8种重金属元素的空间分布规律、污染成因分析以及污染源定位方法。资源包内仅含1个docx文档大小2.94MB包含三维地形图、元素丰度等值线图、地形等高线图等配套图件并附有用于绘图和求解三元二次方程组的Matlab源代码。已有215人学习下载适合数学建模竞赛参赛者、环境数据爱好者以及希望提升Matlab绘图与建模能力的读者参考。通过这份文档读者既可以研读一套结构完整的建模思路也可以参考和复用其中代码来绘制元素分布图并理解污染源识别中的微分方程模型。1. 城市表层土壤重金属污染分析从题目数据到Matlab图件交付数学建模A题城市表层土壤重金属污染分析这道题表面上考的是环境评价与统计建模实际淘汰率最高的环节发生在数据读取和出图阶段。题目给出的Excel里通常包含采样点经纬度、功能区分组和As、Cd、Cr、Cu、Hg、Ni、Pb、Zn八种重金属浓度要求完成污染评价、空间分布推断和污染源识别最后连同所有图件一起提交Matlab源代码。很多队伍把时间耗在公式推导上最后用残存的半天时间仓促画图图件风格不统一、插值边界全是空白、代码里硬编码路径评审体验很差。这篇文章按“数据准备 → 评价计算 → 空间插值 → 聚类源解析 → 图件打包”这条完整路径展开给出一套可以直接套用的Matlab实现。2. 数据清洗与土壤重金属污染评价的Matlab实现2.1 用readtable读取采样坐标、功能区标签与八种重金属浓度常见做法是用readtable替代老旧的xlsread。xlsread在处理多工作表、空白单元格和中文列名时经常返回元胞数组后续还要手动转换类型readtable则会自动识别每列的数据类型把经纬度解析为double把功能区解析为categorical或string省掉大量体力活。% 使用detectImportOptions保留Excel原始列名 opts detectImportOptions(采样点数据.xlsx); opts.PreserveVariableNames true; T readtable(采样点数据.xlsx, opts); % 提取经纬度坐标列 x T.经度; y T.纬度; % 八种重金属的列名后续循环画图、计算都要复用 metals {As, Cd, Cr, Cu, Hg, Ni, Pb, Zn}; Z zeros(height(T), length(metals)); for i 1:length(metals) Z(:, i) T.(metals{i}); end % 功能区分组标签污染评价后按组对比时使用 zone T.功能区;PreserveVariableNames是这里最值得留意的参数。Excel里的中文列名如果不加这个选项读取时可能被替换成Var1这类占位符后面T.(metals{i})的动态列索引就会失败。采样点数据的高度不应该写死用height(T)动态获取这样即使Excel里删掉了几行异常记录代码也不会崩。metals这个元胞数组从始至终只维护一次后续的评价计算、批量画图和聚类分析全部由它驱动避免在多个脚本里各自维护一份金属名称列表导致不一致。清洗阶段重点检查两类问题。第一是缺失值直接用isnan定位如果采样量够大就删除对应行如果删除会影响某功能区的样本数用fillmissing按列均值填充。第二是离群值散点图里明显脱离采样区域、浓度高出周围两个数量级的点通常录入错误这类点应当在插值前剔除否则griddata插值结果会在局部形成一个醒目的“火山口”。2.2 单因子污染指数与内梅罗综合指数的函数封装污染评价有两层计算。单因子指数Pi Ci/Si其中Ci是实测浓度Si是评价标准值。标准值的选取直接影响结论常见的做法是采用《土壤环境质量标准》中的筛选值也可以取当地土壤背景值。筛选值偏宽松背景值偏严格选择哪一个要在报告里明确交代。内梅罗综合指数则是对单因子指数的二次综合公式为P_n sqrt((mean(P_i)^2 max(P_i)^2) / 2)它的设计意图是让最大污染因子在综合评价中占据话语权避免某个重金属严重超标却因平均值较低而被掩盖。这个特性对城市土壤污染场景非常合适因为城市土壤往往是点源污染某一种重金属的局部高值恰恰是需要优先关注的信号。function [Pi, Pn] nemeo_eval(C, S) % C: n行m列实测浓度矩阵行为采样点列为重金属 % S: 1行m列标准值与metals顺序一致 Pi C ./ S; Pm mean(Pi, 2); Pmax max(Pi, [], 2); Pn sqrt((Pm.^2 Pmax.^2) / 2); end这里有两个初学者容易写错的细节。max(Pi)在Matlab里默认沿第一维求解返回的是每个列的最大值需要对每一行求解时必须写成max(Pi, [], 2)第二个维度参数是2第三位缺省时不会做列压缩。mean(Pi, 2)同理mean(Pi)得到的是每个金属的平均值加上2之后才是每个采样点的平均单因子指数。如果需要给不同重金属赋予不同权重可以在Pm处改为加权平均用层次分析法构造判断矩阵计算特征向量后归一化得到权重向量w然后Pm Pi * w(:)这就是在线权重计算中常见的层次分析法matlab代码实现。提示评价标准向量S不要直接写在函数里应当作为参数传入。不同城市、不同题目的评价标准差异很大硬编码会让代码复用性变差。2.3 用功能区作分组统计并输出评价结果表题目数据里一般带有功能区列通常划分为生活区、工业区、主干道、公园绿地和山区。计算完内梅罗指数后按功能区统计均值可以得到“工业区 主干道 生活区 公园绿地 山区”这类评阅专家熟悉的结论链条。[G, zoneName] findgroups(zone); Pn_mean splitapply(mean, Pn, G); Pn_max splitapply(max, Pn, G); outTable table(zoneName, Pn_mean, Pn_max); writetable(outTable, 内梅罗指数_分区统计.xlsx);findgroups返回两个输出第一个是每个样本所属组的编号第二个是这个功能区的中文名称后者直接作为表格行标签很方便。splitapply把Pn按G分组后分别调用mean和max不需要手写for循环。如果还要输出单因子指数的分区均值把Pi用同样的方式传入即可。内梅罗指数的分级标准不同文献略有差异常用的分级如下等级内梅罗指数污染程度IPn ≤ 0.7无污染II0.7 Pn ≤ 1.0尚清洁III1.0 Pn ≤ 2.0轻度污染IV2.0 Pn ≤ 3.0中度污染VPn 3.0重度污染分级阈值在论文里要标注来源。这里0.7和1.0的取值来自国内土壤环境质量评价常用标准部分文献以0.7作为警戒线也有直接以1.0为起点的简化做法评审更看重你是否交代了取值的依据。将分级表与分区统计结果合并输出一张包含“功能区、平均内梅罗指数、污染等级”的汇总表就能作为论文中的评价结果表。3. 用Matlab画图与空间插值生成重金属浓度分布图件3.1 先画散点图排查异常值与采样布局拿到坐标和浓度数据后不要急于插值先画一张散点图。散点图能同时暴露两类问题一是采样点是否存在明显聚集或空洞二是浓度值是否存在录入错误造成的异常高值。figure; scatter(x, y, 40, Z(:, 7), filled); % 40为散点直径filled表示实心 colorbar; colormap(parula(256)); title(Pb浓度空间散点分布); xlabel(经度); ylabel(纬度);scatter的第三个参数是点的直径单位是磅与坐标轴范围无关屏幕尺寸不同时视觉效果会有差别第四个参数是颜色数据每个采样点根据Pb浓度映射到colormap上。parula是Matlab默认的颜色映射对色觉障碍者友好且不会有jet那种从蓝到红的视觉跳跃。如果某个点的颜色明显深于周边且浓度高出两个数量级先回到Excel核对原始记录不要直接带入插值。再配合直方图观察浓度分布形态histogram(Z(:, 7), 30); xlabel(Pb浓度/(mg/kg)); ylabel(频数);如果直方图呈现明显的长尾分布说明数据不符合正态假设后续做主成分分析或聚类前需要取对数变换。城市土壤重金属数据几乎都有这个特征点位靠近污染源的区域浓度成倍偏高直方图右拖尾严重。对数变换比zscore标准化更适合这类偏态数据。3.2 用griddata把离散采样点插值成规则网格散点图只能表达采样点位置的浓度无法支撑“整个城区污染分布”这种面状结论必须做空间插值。Matlab里最直接的插值工具是griddata它能将不规则分布的散点插值到规则网格上。% 生成规则网格200×150是经纬度方向的分辨率 [xg, yg] meshgrid(linspace(min(x), max(x), 200), ... linspace(min(y), max(y), 150)); % v4是matlab画图里比较稳的插值算法边界不易出现大范围NaN zg griddata(x, y, Z(:, 7), xg, yg, v4); % 插值结果边界偶尔会出现NaN用全图最小值兜底 zg(isnan(zg)) min(zg(:));griddata的方法参数需要按数据情况选择。linear速度最快但凸包以外的区域全部返回NaN当采样点在城区边缘呈凹形分布时凹进去的那块多边形会变成白洞natural基于自然邻域插值在Matlab较新版本中可用但数据量大时耗时明显v4基于薄板样条对采样点数量不多、分布不规则的城市土壤数据最稳。网格密度200×150对应四万个插值点对一张图来说足够平滑再加密到400×300对视觉提升有限只会让print导出时文件变大、计算变慢。如果竞赛要求更高可以用fitrgp做高斯过程回归插值相当于给浓度场附加一个不确定性估计在图件上叠加置信区间。作为备选方案griddata的速度优势明显且代码量最少。3.3 用contourf与imagesc输出二维浓度分布图插值完成后绘制填充等值线图。这一步是整张图件最核心的环节很多队伍的图看起来“业余”问题出在没有设置LineStyle和坐标轴方向。figure(Color, w); contourf(xg, yg, zg, 24, LineStyle, none); set(gca, YDir, normal); hold on; plot(x, y, k., MarkerSize, 4); % 叠加原始采样点 colorbar; axis equal tight; title(Pb浓度空间分布);contourf的第四个参数24表示等值线层级数层级越多颜色过渡越细腻但过多的分层会让色差区分度下降一般20到30层足够。LineStyle, none这个参数很容易被忽略默认情况下contourf会在等值线之间画细线导出的图在低分辨率下会出现密集的黑线网格视觉上非常脏。set(gca, YDir, normal)是另一个高频坑点。contourf默认的Y轴方向是反转的也就是纬度小的地方在图像上方这跟地理直觉相反加上这行代码让纬度方向朝上。axis equal保证经度和纬度方向上的物理长度比例一致tight将坐标范围裁剪到插值网格覆盖的区域避免四个方向出现大片空白。叠加采样点用plot(x, y, k.)黑色小点既能表达原始数据的分布密度又不至于遮盖浓度场的色彩信息。作为对比imagesc也能画出类似的热图但它把数据矩阵当作图像像素输出坐标轴变成行列索引需要手动set(gca, XTick, ...)映射经纬度而且高宽比很难控制一般不建议在地理空间分布图里使用。3.4 批量绘制八种重金属图件并按规则导出题目要求提交所有图件意味着八种重金属每样一张浓度分布图加上评价等级图、聚类图总数超过十张。逐张手写脚本既慢又容易漏正确做法是循环驱动。mkdir(figs); for i 1:length(metals) % 对当前金属做插值 zg griddata(x, y, Z(:, i), xg, yg, v4); zg(isnan(zg)) min(zg(:)); % 关闭窗口显示后台渲染大幅加快循环 figure(Visible, off); contourf(xg, yg, zg, 24, LineStyle, none); set(gca, YDir, normal); colorbar; title([metals{i} 浓度空间分布]); axis equal tight; % 导出300dpi PNG提交给评审的图件一般要求清晰 print(gcf, sprintf(figs/%s_浓度空间分布.png, metals{i}), -dpng, -r300); close(gcf); endfigure(Visible, off)让窗口在后台渲染不弹出GUI界面循环跑几十张图时速度提升明显。每次循环结束用close(gcf)释放当前图形句柄否则内存会累积到让Matlab卡顿。文件名采用金属名_浓度空间分布.png这种规则与论文中的图标题一一对应评审对照图表时不需要猜测内容。-r300指定300dpi分辨率电脑屏幕看图没有差别但打印纸质论文时150dpi以上的差异是肉眼可辨的。4. 聚类与相关性分析kmeans与层次聚类识别污染源4.1 标准化与相关系数矩阵先看元素是否共生不同重金属的浓度量纲相同但数值范围差异巨大Hg的浓度可能是零点几毫克每千克而Zn可以到数百如果不做标准化直接算距离Hg的贡献会被完全淹没。常见的做法是zscore标准化将每列变换为零均值、单位方差。Zn zscore(Z); % 相关系数矩阵可视化重金属之间的共生关系 R corr(Zn, Rows, pairwise); figure; imagesc(R); colorbar; axis square; set(gca, XTick, 1:8, XTickLabel, metals); set(gca, YTick, 1:8, YTickLabel, metals);corr计算Pearson相关系数Rows, pairwise表示两个变量做相关性计算时只使用双方都不缺失的样本保留更多有效数据。相关系数矩阵热图能直观反映哪些重金属倾向于共同升高比如工业区常见的Cu、Zn、Pb组合如果相关系数超过0.7意味着它们可能来自同一类污染源。在相关性分析之后紧接着做主成分分析可以将八个变量压缩为两到三个主成分与聚类结果互相验证[coeff, score, latent] pca(Zn); explained latent ./ sum(latent) * 100; % 观察前两个主成分的贡献率及载荷 coeff(:, 1:3)coeff中每个主成分列向量的绝对值大小反映了该主成分主要受哪些重金属驱动latent是特征值。一般情况下前两个主成分的累计贡献率达到60%以上就可以绘制载荷散点图横纵坐标分别是第一、二主成分的载荷对角线方向分布的重金属属于同一污染源。这一步属于matlab图像处理基本操作但它在论文里的说服力很强评阅专家看到PCA载荷图和聚类结果一致时对源解析结论的信任度会明显提高。4.2 kmeans聚类算法matlab对采样点聚类并确定k值在污染源识别中kmeans聚类的典型用法是将采样点按重金属浓度特征分成若干组每组在空间上聚集在一起时对应一个污染源。聚类前必须使用zscore标准化后的数据直接对原始浓度做kmeans时高浓度金属主导聚类结果失去多金属综合语义。rng(42); % 固定随机种子保证结果可复现 eva evalclusters(Zn, kmeans, silhouette, KList, 2:8); k eva.OptimalK; disp([轮廓系数法推荐k num2str(k)]); [idx, C] kmeans(Zn, k, Replicates, 20, MaxIter, 500);evalclusters是确定类别数的标准工具silhouette方法计算每个样本的轮廓系数最优k对应的平均轮廓系数最大。Replicates, 20是kmeans聚类算法matlab实现中最重要的参数kmeans的初值随机选择容易陷入局部最优重复20次后取目标函数最小的结果稳定性显著提升。MaxIter, 500限制每次迭代的次数防止数据量大时单次运行时间过长。聚类结束后用gscatter把聚类标签绘制到经纬度图上figure; gscatter(x, y, idx, [], [], 20); xlabel(经度); ylabel(纬度); title(kmeans聚类结果空间分布); legend(类别1, 类别2, 类别3);如果聚类结果在空间上呈现明显的团块结构比如一类集中在工业区周边、另一类沿主干道延伸就可以结合功能区信息解释污染来源。常见误用是做了聚类却不看空间分布单纯在表格里报一个轮廓系数评阅专家无法判断聚类结果是否有实际地理意义。4.3 用层次聚类与变量聚类热图相关距离怎么设kmeans处理的是采样点对变量也就是八种重金属做层次聚类能得到元素之间的亲疏关系树状图与PCA载荷图互相印证。D pdist(Zn, correlation); L linkage(D, average); figure; dendrogram(L, 0);pdist计算样本间的两两距离correlation将两个样本看作向量后计算它们的相关系数转换为距离。这样选择的原因如果两种重金属在空间上同升同降它们对距离的贡献就小倾向于聚到同一分支这是污染源共生关系最直接的表达。linkage中的average表示类间距离取两类样本间所有距离的均值对异常值稳健ward则以组内方差最小化为目标聚类结果更紧凑适合污染样本这种近似连续过渡的数据。dendrogram(L, 0)的第二个参数为0表示显示全部叶节点八种金属的树状图足够清晰不需要截断。如果需要将层次聚类与热图结合直接调用clustergramclustergram(Zn, ColumnLabels, metals);clustergram会同时显示样本和变量的聚类树HeatMap的行列自动按聚类结果重排相关性高的金属相邻排列视觉上比单纯的热图更有结构感。变量间的距离度量同样推荐correlation如果发现Hg单独一支、Cr与Ni靠近、Pb与Zn与Cu绑定这种结构就可以写入污染来源分析。4.4 用BP神经网络拟合曲线预测未知点浓度作为评价和聚类之外的扩展模型BP神经网络经常被用来做未知点的浓度预测。它的思路是拿经纬度作输入、某重金属浓度为输出训练一个前馈网络再用训练好的网络对规则网格上的待测点做预测效果等同于一种非线性插值。这种做法对竞赛、论文是加分项但要注意它不具备污染源解析的解释能力。X [x, y]; % 输入特征经纬度 Y Z(:, 5); % 以Hg为例 net fitnet([12 8]); % 两层隐层神经元12和8 [net, tr] train(net, X, Y); % 在原始数据上的拟合值与实测对比 yhat net(X); figure; plot(1:length(Y), Y, o, 1:length(Y), yhat, *-); legend(实测, 模型预测);fitnet的第一参数是隐层结构向量[12 8]表示两个隐层神经元数分别是12和8。隐层规模越小拟合能力越弱越大越容易过拟合这个规模对几十个样本量级的数据够用。训练函数默认是trainlm即Levenberg-Marquardt收敛速度快但占用内存高数据集超过几百个样本时建议切换为trainscg。BP神经网络拟合曲线适合用来论证“未知坐标点浓度可预测”这个结论但脱离经纬度之外如果有海拔、距工厂距离等协变量应一并加入输入特征否则纯靠经纬度插值的意义与griddata差别不大。5. 图件统一风格与Matlab源代码交付的打包技巧5.1 用matlab画图统一字体与dpi的通用配置评审拿到十几张图件时最先感受到的是整体观感。字号忽大忽小、背景色一白一灰、某些图的色标刻度线缺失都会拉低成品的专业度。推荐在出图脚本开头一次性设置全局默认值set(0, DefaultAxesFontName, 宋体); set(0, DefaultAxesFontSize, 11); set(0, DefaultFigureColor, w); set(0, DefaultAxesLineWidth, 1.2); set(0, DefaultTextInterpreter, none);set(0, ...)设置的是root级别的默认属性对之后创建的所有figure生效。其中DefaultTextInterpreter最容易被忽略默认值是tex元素名称中的下划线会被解释成下标比如Cu_浓度会因为下划线显示异常改成none后文本按字面显示。字号11对应Word正文的小四打印出来比例协调。DefaultFigureColor设为w保证导出的PNG背景为纯白避免灰色背景在Word页面里显得突兀。图件导出规格建议按下表执行用途格式分辨率说明论文插图PNG300dpi体积可控插入Word清晰打印答辩TIFF600dpi文件大一般最后提交时转换二次修改FIG矢量保留图层与数据方便调整在线预览PDF矢量合并多图输出便于评委翻页5.2 从变量名自动生成标题与文件名整理可运行的交付包图件批量化输出的最后一步是用变量自动拼接标题和文件名避免手写八个title和八个saveas。配合exportgraphics还能把多张图合并为一个PDF便于评委连续翻阅regionName 主城区; for i 1:length(metals) zg griddata(x, y, Z(:, i), xg, yg, v4); figure(Visible, off); contourf(xg, yg, zg, 24, LineStyle, none); colorbar; title(sprintf(%s土壤%s浓度空间分布, regionName, metals{i})); axis equal tight; exportgraphics(gcf, sprintf(图件/%s_%s分布.pdf, regionName, metals{i}), Resolution, 300); close(gcf); endexportgraphics是比print更现代的导出方式它的特点是导出的范围精确对应坐标区不会像print那样把Figure的空白边也打进去多条矢量字体保留完整支持Append, true参数实现多页PDF追加。print在导出单张PNG时仍然够用但涉及多页报告或对边缘裁剪有要求时exportgraphics更合适。源代码交付时至少拆成四个脚本read_data.m只负责读Excel和清洗evaluation.m计算单因子与内梅罗指数并输出表格interpolation.m完成网格插值plot_all.m统一输出所有图件。每个脚本头部用两到三行注释写明依赖的前序脚本和数据文件名不要出现本地绝对路径使用which(采样点数据.xlsx)定位文件所在目录。运行plot_all.m能重新生成全部图件并且与论文里的图号一一对应这份交付物才算可用。本文还有配套的精品资源点击获取
Matlab实现城市土壤重金属污染分析与空间插值出图实战
简介一份完整的数学建模竞赛A题参赛文档以城市表层土壤重金属污染分析为项目背景系统讲解基于Matlab的环境数据建模与可视化流程。内容涵盖题目重述、问题分析、模型假设、模型建立与求解、模型评价与推广等完整环节重点展示As、Cd、Cr、Cu、Hg、Ni、Pb、Zn这8种重金属元素的空间分布规律、污染成因分析以及污染源定位方法。资源包内仅含1个docx文档大小2.94MB包含三维地形图、元素丰度等值线图、地形等高线图等配套图件并附有用于绘图和求解三元二次方程组的Matlab源代码。已有215人学习下载适合数学建模竞赛参赛者、环境数据爱好者以及希望提升Matlab绘图与建模能力的读者参考。通过这份文档读者既可以研读一套结构完整的建模思路也可以参考和复用其中代码来绘制元素分布图并理解污染源识别中的微分方程模型。1. 城市表层土壤重金属污染分析从题目数据到Matlab图件交付数学建模A题城市表层土壤重金属污染分析这道题表面上考的是环境评价与统计建模实际淘汰率最高的环节发生在数据读取和出图阶段。题目给出的Excel里通常包含采样点经纬度、功能区分组和As、Cd、Cr、Cu、Hg、Ni、Pb、Zn八种重金属浓度要求完成污染评价、空间分布推断和污染源识别最后连同所有图件一起提交Matlab源代码。很多队伍把时间耗在公式推导上最后用残存的半天时间仓促画图图件风格不统一、插值边界全是空白、代码里硬编码路径评审体验很差。这篇文章按“数据准备 → 评价计算 → 空间插值 → 聚类源解析 → 图件打包”这条完整路径展开给出一套可以直接套用的Matlab实现。2. 数据清洗与土壤重金属污染评价的Matlab实现2.1 用readtable读取采样坐标、功能区标签与八种重金属浓度常见做法是用readtable替代老旧的xlsread。xlsread在处理多工作表、空白单元格和中文列名时经常返回元胞数组后续还要手动转换类型readtable则会自动识别每列的数据类型把经纬度解析为double把功能区解析为categorical或string省掉大量体力活。% 使用detectImportOptions保留Excel原始列名 opts detectImportOptions(采样点数据.xlsx); opts.PreserveVariableNames true; T readtable(采样点数据.xlsx, opts); % 提取经纬度坐标列 x T.经度; y T.纬度; % 八种重金属的列名后续循环画图、计算都要复用 metals {As, Cd, Cr, Cu, Hg, Ni, Pb, Zn}; Z zeros(height(T), length(metals)); for i 1:length(metals) Z(:, i) T.(metals{i}); end % 功能区分组标签污染评价后按组对比时使用 zone T.功能区;PreserveVariableNames是这里最值得留意的参数。Excel里的中文列名如果不加这个选项读取时可能被替换成Var1这类占位符后面T.(metals{i})的动态列索引就会失败。采样点数据的高度不应该写死用height(T)动态获取这样即使Excel里删掉了几行异常记录代码也不会崩。metals这个元胞数组从始至终只维护一次后续的评价计算、批量画图和聚类分析全部由它驱动避免在多个脚本里各自维护一份金属名称列表导致不一致。清洗阶段重点检查两类问题。第一是缺失值直接用isnan定位如果采样量够大就删除对应行如果删除会影响某功能区的样本数用fillmissing按列均值填充。第二是离群值散点图里明显脱离采样区域、浓度高出周围两个数量级的点通常录入错误这类点应当在插值前剔除否则griddata插值结果会在局部形成一个醒目的“火山口”。2.2 单因子污染指数与内梅罗综合指数的函数封装污染评价有两层计算。单因子指数Pi Ci/Si其中Ci是实测浓度Si是评价标准值。标准值的选取直接影响结论常见的做法是采用《土壤环境质量标准》中的筛选值也可以取当地土壤背景值。筛选值偏宽松背景值偏严格选择哪一个要在报告里明确交代。内梅罗综合指数则是对单因子指数的二次综合公式为P_n sqrt((mean(P_i)^2 max(P_i)^2) / 2)它的设计意图是让最大污染因子在综合评价中占据话语权避免某个重金属严重超标却因平均值较低而被掩盖。这个特性对城市土壤污染场景非常合适因为城市土壤往往是点源污染某一种重金属的局部高值恰恰是需要优先关注的信号。function [Pi, Pn] nemeo_eval(C, S) % C: n行m列实测浓度矩阵行为采样点列为重金属 % S: 1行m列标准值与metals顺序一致 Pi C ./ S; Pm mean(Pi, 2); Pmax max(Pi, [], 2); Pn sqrt((Pm.^2 Pmax.^2) / 2); end这里有两个初学者容易写错的细节。max(Pi)在Matlab里默认沿第一维求解返回的是每个列的最大值需要对每一行求解时必须写成max(Pi, [], 2)第二个维度参数是2第三位缺省时不会做列压缩。mean(Pi, 2)同理mean(Pi)得到的是每个金属的平均值加上2之后才是每个采样点的平均单因子指数。如果需要给不同重金属赋予不同权重可以在Pm处改为加权平均用层次分析法构造判断矩阵计算特征向量后归一化得到权重向量w然后Pm Pi * w(:)这就是在线权重计算中常见的层次分析法matlab代码实现。提示评价标准向量S不要直接写在函数里应当作为参数传入。不同城市、不同题目的评价标准差异很大硬编码会让代码复用性变差。2.3 用功能区作分组统计并输出评价结果表题目数据里一般带有功能区列通常划分为生活区、工业区、主干道、公园绿地和山区。计算完内梅罗指数后按功能区统计均值可以得到“工业区 主干道 生活区 公园绿地 山区”这类评阅专家熟悉的结论链条。[G, zoneName] findgroups(zone); Pn_mean splitapply(mean, Pn, G); Pn_max splitapply(max, Pn, G); outTable table(zoneName, Pn_mean, Pn_max); writetable(outTable, 内梅罗指数_分区统计.xlsx);findgroups返回两个输出第一个是每个样本所属组的编号第二个是这个功能区的中文名称后者直接作为表格行标签很方便。splitapply把Pn按G分组后分别调用mean和max不需要手写for循环。如果还要输出单因子指数的分区均值把Pi用同样的方式传入即可。内梅罗指数的分级标准不同文献略有差异常用的分级如下等级内梅罗指数污染程度IPn ≤ 0.7无污染II0.7 Pn ≤ 1.0尚清洁III1.0 Pn ≤ 2.0轻度污染IV2.0 Pn ≤ 3.0中度污染VPn 3.0重度污染分级阈值在论文里要标注来源。这里0.7和1.0的取值来自国内土壤环境质量评价常用标准部分文献以0.7作为警戒线也有直接以1.0为起点的简化做法评审更看重你是否交代了取值的依据。将分级表与分区统计结果合并输出一张包含“功能区、平均内梅罗指数、污染等级”的汇总表就能作为论文中的评价结果表。3. 用Matlab画图与空间插值生成重金属浓度分布图件3.1 先画散点图排查异常值与采样布局拿到坐标和浓度数据后不要急于插值先画一张散点图。散点图能同时暴露两类问题一是采样点是否存在明显聚集或空洞二是浓度值是否存在录入错误造成的异常高值。figure; scatter(x, y, 40, Z(:, 7), filled); % 40为散点直径filled表示实心 colorbar; colormap(parula(256)); title(Pb浓度空间散点分布); xlabel(经度); ylabel(纬度);scatter的第三个参数是点的直径单位是磅与坐标轴范围无关屏幕尺寸不同时视觉效果会有差别第四个参数是颜色数据每个采样点根据Pb浓度映射到colormap上。parula是Matlab默认的颜色映射对色觉障碍者友好且不会有jet那种从蓝到红的视觉跳跃。如果某个点的颜色明显深于周边且浓度高出两个数量级先回到Excel核对原始记录不要直接带入插值。再配合直方图观察浓度分布形态histogram(Z(:, 7), 30); xlabel(Pb浓度/(mg/kg)); ylabel(频数);如果直方图呈现明显的长尾分布说明数据不符合正态假设后续做主成分分析或聚类前需要取对数变换。城市土壤重金属数据几乎都有这个特征点位靠近污染源的区域浓度成倍偏高直方图右拖尾严重。对数变换比zscore标准化更适合这类偏态数据。3.2 用griddata把离散采样点插值成规则网格散点图只能表达采样点位置的浓度无法支撑“整个城区污染分布”这种面状结论必须做空间插值。Matlab里最直接的插值工具是griddata它能将不规则分布的散点插值到规则网格上。% 生成规则网格200×150是经纬度方向的分辨率 [xg, yg] meshgrid(linspace(min(x), max(x), 200), ... linspace(min(y), max(y), 150)); % v4是matlab画图里比较稳的插值算法边界不易出现大范围NaN zg griddata(x, y, Z(:, 7), xg, yg, v4); % 插值结果边界偶尔会出现NaN用全图最小值兜底 zg(isnan(zg)) min(zg(:));griddata的方法参数需要按数据情况选择。linear速度最快但凸包以外的区域全部返回NaN当采样点在城区边缘呈凹形分布时凹进去的那块多边形会变成白洞natural基于自然邻域插值在Matlab较新版本中可用但数据量大时耗时明显v4基于薄板样条对采样点数量不多、分布不规则的城市土壤数据最稳。网格密度200×150对应四万个插值点对一张图来说足够平滑再加密到400×300对视觉提升有限只会让print导出时文件变大、计算变慢。如果竞赛要求更高可以用fitrgp做高斯过程回归插值相当于给浓度场附加一个不确定性估计在图件上叠加置信区间。作为备选方案griddata的速度优势明显且代码量最少。3.3 用contourf与imagesc输出二维浓度分布图插值完成后绘制填充等值线图。这一步是整张图件最核心的环节很多队伍的图看起来“业余”问题出在没有设置LineStyle和坐标轴方向。figure(Color, w); contourf(xg, yg, zg, 24, LineStyle, none); set(gca, YDir, normal); hold on; plot(x, y, k., MarkerSize, 4); % 叠加原始采样点 colorbar; axis equal tight; title(Pb浓度空间分布);contourf的第四个参数24表示等值线层级数层级越多颜色过渡越细腻但过多的分层会让色差区分度下降一般20到30层足够。LineStyle, none这个参数很容易被忽略默认情况下contourf会在等值线之间画细线导出的图在低分辨率下会出现密集的黑线网格视觉上非常脏。set(gca, YDir, normal)是另一个高频坑点。contourf默认的Y轴方向是反转的也就是纬度小的地方在图像上方这跟地理直觉相反加上这行代码让纬度方向朝上。axis equal保证经度和纬度方向上的物理长度比例一致tight将坐标范围裁剪到插值网格覆盖的区域避免四个方向出现大片空白。叠加采样点用plot(x, y, k.)黑色小点既能表达原始数据的分布密度又不至于遮盖浓度场的色彩信息。作为对比imagesc也能画出类似的热图但它把数据矩阵当作图像像素输出坐标轴变成行列索引需要手动set(gca, XTick, ...)映射经纬度而且高宽比很难控制一般不建议在地理空间分布图里使用。3.4 批量绘制八种重金属图件并按规则导出题目要求提交所有图件意味着八种重金属每样一张浓度分布图加上评价等级图、聚类图总数超过十张。逐张手写脚本既慢又容易漏正确做法是循环驱动。mkdir(figs); for i 1:length(metals) % 对当前金属做插值 zg griddata(x, y, Z(:, i), xg, yg, v4); zg(isnan(zg)) min(zg(:)); % 关闭窗口显示后台渲染大幅加快循环 figure(Visible, off); contourf(xg, yg, zg, 24, LineStyle, none); set(gca, YDir, normal); colorbar; title([metals{i} 浓度空间分布]); axis equal tight; % 导出300dpi PNG提交给评审的图件一般要求清晰 print(gcf, sprintf(figs/%s_浓度空间分布.png, metals{i}), -dpng, -r300); close(gcf); endfigure(Visible, off)让窗口在后台渲染不弹出GUI界面循环跑几十张图时速度提升明显。每次循环结束用close(gcf)释放当前图形句柄否则内存会累积到让Matlab卡顿。文件名采用金属名_浓度空间分布.png这种规则与论文中的图标题一一对应评审对照图表时不需要猜测内容。-r300指定300dpi分辨率电脑屏幕看图没有差别但打印纸质论文时150dpi以上的差异是肉眼可辨的。4. 聚类与相关性分析kmeans与层次聚类识别污染源4.1 标准化与相关系数矩阵先看元素是否共生不同重金属的浓度量纲相同但数值范围差异巨大Hg的浓度可能是零点几毫克每千克而Zn可以到数百如果不做标准化直接算距离Hg的贡献会被完全淹没。常见的做法是zscore标准化将每列变换为零均值、单位方差。Zn zscore(Z); % 相关系数矩阵可视化重金属之间的共生关系 R corr(Zn, Rows, pairwise); figure; imagesc(R); colorbar; axis square; set(gca, XTick, 1:8, XTickLabel, metals); set(gca, YTick, 1:8, YTickLabel, metals);corr计算Pearson相关系数Rows, pairwise表示两个变量做相关性计算时只使用双方都不缺失的样本保留更多有效数据。相关系数矩阵热图能直观反映哪些重金属倾向于共同升高比如工业区常见的Cu、Zn、Pb组合如果相关系数超过0.7意味着它们可能来自同一类污染源。在相关性分析之后紧接着做主成分分析可以将八个变量压缩为两到三个主成分与聚类结果互相验证[coeff, score, latent] pca(Zn); explained latent ./ sum(latent) * 100; % 观察前两个主成分的贡献率及载荷 coeff(:, 1:3)coeff中每个主成分列向量的绝对值大小反映了该主成分主要受哪些重金属驱动latent是特征值。一般情况下前两个主成分的累计贡献率达到60%以上就可以绘制载荷散点图横纵坐标分别是第一、二主成分的载荷对角线方向分布的重金属属于同一污染源。这一步属于matlab图像处理基本操作但它在论文里的说服力很强评阅专家看到PCA载荷图和聚类结果一致时对源解析结论的信任度会明显提高。4.2 kmeans聚类算法matlab对采样点聚类并确定k值在污染源识别中kmeans聚类的典型用法是将采样点按重金属浓度特征分成若干组每组在空间上聚集在一起时对应一个污染源。聚类前必须使用zscore标准化后的数据直接对原始浓度做kmeans时高浓度金属主导聚类结果失去多金属综合语义。rng(42); % 固定随机种子保证结果可复现 eva evalclusters(Zn, kmeans, silhouette, KList, 2:8); k eva.OptimalK; disp([轮廓系数法推荐k num2str(k)]); [idx, C] kmeans(Zn, k, Replicates, 20, MaxIter, 500);evalclusters是确定类别数的标准工具silhouette方法计算每个样本的轮廓系数最优k对应的平均轮廓系数最大。Replicates, 20是kmeans聚类算法matlab实现中最重要的参数kmeans的初值随机选择容易陷入局部最优重复20次后取目标函数最小的结果稳定性显著提升。MaxIter, 500限制每次迭代的次数防止数据量大时单次运行时间过长。聚类结束后用gscatter把聚类标签绘制到经纬度图上figure; gscatter(x, y, idx, [], [], 20); xlabel(经度); ylabel(纬度); title(kmeans聚类结果空间分布); legend(类别1, 类别2, 类别3);如果聚类结果在空间上呈现明显的团块结构比如一类集中在工业区周边、另一类沿主干道延伸就可以结合功能区信息解释污染来源。常见误用是做了聚类却不看空间分布单纯在表格里报一个轮廓系数评阅专家无法判断聚类结果是否有实际地理意义。4.3 用层次聚类与变量聚类热图相关距离怎么设kmeans处理的是采样点对变量也就是八种重金属做层次聚类能得到元素之间的亲疏关系树状图与PCA载荷图互相印证。D pdist(Zn, correlation); L linkage(D, average); figure; dendrogram(L, 0);pdist计算样本间的两两距离correlation将两个样本看作向量后计算它们的相关系数转换为距离。这样选择的原因如果两种重金属在空间上同升同降它们对距离的贡献就小倾向于聚到同一分支这是污染源共生关系最直接的表达。linkage中的average表示类间距离取两类样本间所有距离的均值对异常值稳健ward则以组内方差最小化为目标聚类结果更紧凑适合污染样本这种近似连续过渡的数据。dendrogram(L, 0)的第二个参数为0表示显示全部叶节点八种金属的树状图足够清晰不需要截断。如果需要将层次聚类与热图结合直接调用clustergramclustergram(Zn, ColumnLabels, metals);clustergram会同时显示样本和变量的聚类树HeatMap的行列自动按聚类结果重排相关性高的金属相邻排列视觉上比单纯的热图更有结构感。变量间的距离度量同样推荐correlation如果发现Hg单独一支、Cr与Ni靠近、Pb与Zn与Cu绑定这种结构就可以写入污染来源分析。4.4 用BP神经网络拟合曲线预测未知点浓度作为评价和聚类之外的扩展模型BP神经网络经常被用来做未知点的浓度预测。它的思路是拿经纬度作输入、某重金属浓度为输出训练一个前馈网络再用训练好的网络对规则网格上的待测点做预测效果等同于一种非线性插值。这种做法对竞赛、论文是加分项但要注意它不具备污染源解析的解释能力。X [x, y]; % 输入特征经纬度 Y Z(:, 5); % 以Hg为例 net fitnet([12 8]); % 两层隐层神经元12和8 [net, tr] train(net, X, Y); % 在原始数据上的拟合值与实测对比 yhat net(X); figure; plot(1:length(Y), Y, o, 1:length(Y), yhat, *-); legend(实测, 模型预测);fitnet的第一参数是隐层结构向量[12 8]表示两个隐层神经元数分别是12和8。隐层规模越小拟合能力越弱越大越容易过拟合这个规模对几十个样本量级的数据够用。训练函数默认是trainlm即Levenberg-Marquardt收敛速度快但占用内存高数据集超过几百个样本时建议切换为trainscg。BP神经网络拟合曲线适合用来论证“未知坐标点浓度可预测”这个结论但脱离经纬度之外如果有海拔、距工厂距离等协变量应一并加入输入特征否则纯靠经纬度插值的意义与griddata差别不大。5. 图件统一风格与Matlab源代码交付的打包技巧5.1 用matlab画图统一字体与dpi的通用配置评审拿到十几张图件时最先感受到的是整体观感。字号忽大忽小、背景色一白一灰、某些图的色标刻度线缺失都会拉低成品的专业度。推荐在出图脚本开头一次性设置全局默认值set(0, DefaultAxesFontName, 宋体); set(0, DefaultAxesFontSize, 11); set(0, DefaultFigureColor, w); set(0, DefaultAxesLineWidth, 1.2); set(0, DefaultTextInterpreter, none);set(0, ...)设置的是root级别的默认属性对之后创建的所有figure生效。其中DefaultTextInterpreter最容易被忽略默认值是tex元素名称中的下划线会被解释成下标比如Cu_浓度会因为下划线显示异常改成none后文本按字面显示。字号11对应Word正文的小四打印出来比例协调。DefaultFigureColor设为w保证导出的PNG背景为纯白避免灰色背景在Word页面里显得突兀。图件导出规格建议按下表执行用途格式分辨率说明论文插图PNG300dpi体积可控插入Word清晰打印答辩TIFF600dpi文件大一般最后提交时转换二次修改FIG矢量保留图层与数据方便调整在线预览PDF矢量合并多图输出便于评委翻页5.2 从变量名自动生成标题与文件名整理可运行的交付包图件批量化输出的最后一步是用变量自动拼接标题和文件名避免手写八个title和八个saveas。配合exportgraphics还能把多张图合并为一个PDF便于评委连续翻阅regionName 主城区; for i 1:length(metals) zg griddata(x, y, Z(:, i), xg, yg, v4); figure(Visible, off); contourf(xg, yg, zg, 24, LineStyle, none); colorbar; title(sprintf(%s土壤%s浓度空间分布, regionName, metals{i})); axis equal tight; exportgraphics(gcf, sprintf(图件/%s_%s分布.pdf, regionName, metals{i}), Resolution, 300); close(gcf); endexportgraphics是比print更现代的导出方式它的特点是导出的范围精确对应坐标区不会像print那样把Figure的空白边也打进去多条矢量字体保留完整支持Append, true参数实现多页PDF追加。print在导出单张PNG时仍然够用但涉及多页报告或对边缘裁剪有要求时exportgraphics更合适。源代码交付时至少拆成四个脚本read_data.m只负责读Excel和清洗evaluation.m计算单因子与内梅罗指数并输出表格interpolation.m完成网格插值plot_all.m统一输出所有图件。每个脚本头部用两到三行注释写明依赖的前序脚本和数据文件名不要出现本地绝对路径使用which(采样点数据.xlsx)定位文件所在目录。运行plot_all.m能重新生成全部图件并且与论文里的图号一一对应这份交付物才算可用。本文还有配套的精品资源点击获取