简介面向数学建模竞赛参赛者及环境数据分析学习者的完整资料包以城市表层土壤重金属污染分析A题为载体呈现了从数据预处理、三维地形绘制到污染源定位的全流程建模方案。资源为单个docx文档包体大小约2.94MB内容包含Matlab源代码、等值线图件生成过程、丰度计算方法、基于微分方程与二次方程组求解的污染源推断模型并附有模型评价与推广分析。已有215人学习下载。通过阅读可清晰掌握As、Hg、Zn、Cd、Cr、Cu等多种重金属空间分布的可视化实现理解如何结合风向、水流等条件量化传播特征同时获得可直接复用的绘图脚本与建模思路适合备赛期间快速参考、对照代码练习或作为课程设计完善蓝本。1. 数学建模A题城市表层土壤重金属污染分析真正要交付的是“图件源代码”如果你的电脑上还躺着《数学建模A题城市表层土壤重金属污染分析附所有图件的Matlab源代码.docx》这个文件大概率是当年没做完或者现在要做类似项目却没有完整技术链路。这道题的交付物不是一篇论文而是两样东西能复现的所有图件以及把这些图件忠实生成出来的Matlab源代码。图件要能看出8种重金属在表层土壤中的空间分布、污染等级、污染来源组合和传播趋势代码要能让别人直接跑通。以下就是我梳理的从xlsx到成品图的完整处理链路按数据清洗、空间插值、污染评价、源解析、批量出图五个环节展开所有脚本基于MATLAB R2018b及以上版本。2. 录入Excel到Matlab的第一步读取、清洗和异常值判定2.1 用readtable把采样表变成可计算的tableA题附件通常会提供一个Excel里面是采样点坐标和重金属浓度。我的习惯是直接用readtable读入不手动xlsread因为readtable返回的table类型能保留列名后续用data.As这种写法可读性好得多。假设文件名为soil_data.xlsx工作表叫采样点列名是X,Y,As,Cd,Cr,Cu,Hg,Ni,Pb,Zn那么代码如下raw readtable(soil_data.xlsx, Sheet, 采样点, ... VariableNamingRule, preserve); % 显示前5行确认列名和类型 head(raw, 5) % 查看每列的最小值、最大值、缺失值个数 summary(raw)readtable会自动把Excel里的数值列识别为double但要注意“未检出”会被写成ND或0.01这会导致整列变成cell。这种情况要单独处理用readcell读原始单元格再根据符号过滤或者直接让数据提供方先把ND替换为检测限的1/2。我在多个项目里碰到过这种问题最省事的是在Excel里先处理如果对方做不到就在Matlab里用字符串替换再str2double但那样效率很低还会引入精度丢失。2.2 用盒图和isoutlier识别“高值”与“错误值”的边界土壤重金属数据的特点是天然存在高值这些高值可能来自真实污染点不能看见盒图外的点就删。常见做法是先对每种元素画盒图观察是否有“物理不可能”的值比如负数或明显超过仪器检出上限的数值。有了table之后画盒图非常快metals {As,Cd,Cr,Cu,Hg,Ni,Pb,Zn}; figure(Color,w); for i 1:8 subplot(2,4,i); boxplot(raw.(metals{i}), Symbol, r.); title(metals{i}, FontSize, 11); end sgtitle(原始土壤重金属浓度盒图, FontSize, 14);盒图外部的点只是离群点不全是错误。处理原则是浓度小于0直接删大于3倍四分位距的记下来核对采样记录不要着急剔除。若题目给出了背景值可以把显著高于背景值但远大于样本分布的单个值视为异常用isoutlier指定中位数法tf isoutlier(raw.As, median); raw(tf, :) % 查看这些行isoutlier的median方法用3倍局部中位数绝对偏差作为阈值比默认的quartiles更稳健。这里建议保留这些点因为污染分析需要极端值只有在确认录入错误时才删除。删除要记录哪些行被删了我习惯在代码里加一行日志fprintf(删除异常行: %s\n, join(string(find(tf)), ,));2.3 重复坐标去重与数据一致性校验坐标重复点是土壤数据常遇到的问题。如果同一个X,Y出现了多次插值时会因为函数输入重复导致scatteredInterpolant报错。解决办法是先用unique按行去重[~, ia, ic] unique([X, Y], rows, stable); X X(ia); Y Y(ia); Z Z(ia);ic里保存的是原始行号到去重行号的映射如果你后面要做“每个采样点对应功能区”的关联可以通过ic把去重前的其他字段也同步过来。比如功能区块可以这样提取funcZone raw.FunctionZone(ia);这样就保证了坐标、浓度、功能区三类信息在同一条记录里不会因为去重而错位。3. 从离散采样点到连续污染面空间插值方案怎么选3.1 scatteredInterpolant的插值方法对比Matlab内置scatteredInterpolant支持三种插值方法nearest、linear、natural。我一般不用nearest因为画出来是块状边界不美观linear沿三角网线性插值会产生棱线对污染分布这种连续自然现象来说过于生硬natural基于自然邻域兼顾光滑和保留局部尖峰是重金属污染图的首选。代码如下F_natural scatteredInterpolant(X, Y, Z, natural, none); [Xi, Yi] meshgrid(linspace(min(X), max(X), 150), linspace(min(Y), max(Y), 150)); Zi_natural F_natural(Xi, Yi);这里第二个参数none表示不进行外推因为我们不需要采样范围之外的预测值。如果你选了linear外推四个角落会出现极大的不合理值必须配合掩膜裁剪。如果不填外推方法默认是linear这时候meshgrid生成的四个角会被赋予离谱的数值所以我总是显式写none。插值方法表面特征对尖峰保留计算速度nearest块状台阶差最快linear三角面棱角好快natural光滑最好中对于A题这种100个采样点以内的数据natural方法完全没有性能压力所以优先选它。3.2 网格边界裁剪别让插值跑到采样范围外用boundary函数生成采样点外轮廓把轮廓外的网格值设为NaN这一步能有效避免“地图四角被插值填满”的尴尬k boundary(X, Y, 0.5); % 收缩系数0.5越小越贴合点集 in inpolygon(Xi, Yi, X(k), Y(k)); Zi_natural(~in) NaN;boundary的第三个参数shrink factor控制了包络的紧致程度。值越接近1轮廓越紧致能把凹进去的区域也包进来值越接近0轮廓越接近凸包。对土壤采样点这种分布不均匀的数据0.5通常效果不错。绘制污染分布图figure(Color,w); contourf(Xi, Yi, Zi_natural, 20, LineWidth, 0.5); hold on; scatter(X, Y, 20, Z, filled, MarkerEdgeColor,k); colormap(jet); colorbar; axis equal; xlabel(X (m)); ylabel(Y (m));3.3 如果想用克里金插值的替代方案Matlab没有内置Kriging函数但自己实现变差函数拟合又太繁琐。一个干净的替代方案是用fitrgp高斯过程回归来近似克里金它返回的是正态分布预测拿均值当克里金插值结果正合适gpr fitrgp([X, Y], Z, KernelFunction, squaredExponential); [Zp, ~, interval] predict(gpr, [Xi(:), Yi(:)]); Zp reshape(Zp, size(Xi));fitrgp的优点是能给出置信区间这对以后做“污染不确定性分析”非常有用。不过它的训练时间随样本量增加很快土壤数据一般几百个点还算ok。如果只是画分布图还是natural更顺手。4. 单因子指数与内梅罗指数污染评价的量化指标计算4.1 按功能区设定背景值C0A题附件一般会给出各功能区的背景值表格如生活区、工业区、山区等。不同功能区背景值不同不能拿同一个标准去计算。我的做法是先将功能区列读入然后通过groupsummary按区统计再分别除以对应的背景值。给定背景值向量C0长度8即可算每个采样点的单因子指数% 假设data是清洗后的table包含metals和FunctionZone C0 [25.0, 0.4, 80, 30, 0.2, 30, 50, 100]; % 示例值实际按附件填写 Pij data{:, metals} ./ C0;注意这里除出来是每个采样点每种元素的单因子指数。如果背景值是分功能区的就先把data按功能区分组再分别除以对应的C0不要用for循环用rowfun或splitapply更好。比如生活区的背景值是C0Life工业区是C0Ind可以这样func data.FunctionZone; C0_all zeros(size(data,1), 8); C0_all(func 生活区, :) repmat(C0Life, sum(func生活区), 1); C0_all(func 工业区, :) repmat(C0Ind, sum(func工业区), 1); Pij data{:, metals} ./ C0_all;4.2 计算单项指数与综合指数并映射污染等级内梅罗综合指数强调最大单因子公式是PN sqrt((Pmax^2 Pave^2) / 2)其中Pmax是同一采样点所有元素单因子指数的最大值Pave是平均值。代码Pmax max(Pij, [], 2); Pave mean(Pij, 2); PN sqrt((Pmax.^2 Pave.^2) / 2);然后给每个采样点打上污染等级。常见分级是PN1清洁1~2轻度2~3中度3重度。也有用5级的这里按题目给出的标准来。可以用discretize函数映射到文本edges [0, 1, 2, 3, Inf]; labels {清洁,轻度污染,中度污染,重度污染}; level discretize(PN, edges, categorical, labels);discretize的边界是左闭右开所以PN1会分入1~2的轻度污染。如果想严格按≤算可以把edges设成[0, 1eps, 2eps, ...]避免因为浮点误差导致分界值落入错误区间。4.3 评价结果的空间化把等级画到图上有了level之后把等级标签作为类别用scatter按类别着色得到功能区污染评价图。同时可以统计各功能区不同污染等级的占比用crosstabtab crosstab(data.FunctionZone, level); bar(tab, stacked); legend(categories(level), Location, northwest); grid on;crosstab返回的是交叉频数配合bar画堆叠柱状图这是论文里“不同功能区污染程度对比”的标准图。要注意crosstab的标签顺序是按类别排序的图例顺序要和条形对应上否则会误导读者。5. 用主成分分析反推污染源从协方差到因子旋转5.1 pca函数的用法与载荷矩阵计算主成分分析是识别重金属来源最常用的方法。用Matlab自带pca函数输入是n个采样点×8种重金属的浓度矩阵先做z-score标准化以消除量纲影响stdMetals zscore(data{:, metals}); [coeff, score, latent, ~, explained] pca(stdMetals);coeff是8×8的特征向量矩阵score是主成分得分latent是特征值explained是每个主成分的方差百分比。要得到因子载荷即变量与主成分的相关系数需要将coeff每列乘以对应特征值的平方根loadings coeff .* sqrt(latent);然后看前两个主成分的载荷figure(Color,w); gscatter(loadings(:,1), loadings(:,2), metals, [], oo, 15); xlabel([PC1 (, num2str(round(explained(1))), %)]); ylabel([PC2 (, num2str(round(explained(2))), %)]); grid on;如果只想保留前3个分量做因子旋转可以这样[rotated, T] rotatefactors(loadings(:,1:3), Method, varimax);rotatefactors返回旋转后的载荷矩阵使得每个变量尽可能只在一个因子上有较大载荷方便解释污染源。旋转后的载荷矩阵可以直接打印出来检查哪些元素在同一个因子上载荷高。5.2 载荷图与元素组合的解释套路解释污染源没有标准答案但有一些常见规律Pb和Zn往往来自交通和工业废气Cd和Hg富集于冶炼和燃煤As常与化工或采矿有关。通过载荷图可以把具有相似载荷方向的元素归为一类认为它们同源。例如如果第一因子上Cr和Ni载荷都高可能说明它们来自同一片工业区。注意主成分分析只告诉我们“统计上相关”不能直接证明因果所以论文里的表述应该是“与…相关的排放源”而不是“就是”。5.3 结合K-means做功能区污染聚类如果还想进一步把采样点分类可以用kmeans对score的前3列聚类直接在地图上圈出几块污染聚集区idx kmeans(score(:,1:3), 4, Replicates, 20); gscatter(X, Y, idx);聚类个数4是经验值可以用evalclusters的Davies-Bouldin或Silhouette指标来定eva evalclusters(score(:,1:3), kmeans, silhouette, KList, 2:6);当eva.OptimalK落在3~5之间时说明污染源类型比较稳定。这个方法与主成分分析互补能在地图上圈出几块污染源聚集区是论文里“污染源分布”的支撑图。6. 批量导出300DPI图件并自动拼接Word报告6.1 用exportgraphics统一导出PNGMatlab从R2020a起主推exportgraphics它比print和saveas更稳定能保留透明背景和字体。我写一个循环把前面所有图一次性导出exportgraphics(gcf, sprintf(figure_%02d.png, figIndex), Resolution, 300);注意exportgraphics默认裁剪到图窗中的内容不会像print -dpng那样保留大量空白。建议在导出前统一设置字号set(gcf, Color,w, DefaultAxesFontSize, 12, DefaultTextFontSize, 12);这样所有图件的字号一致报告里不会出现大小不一的丑图。如果你要导出的是contourf记得在导出前调用colormap(jet)否则可能沿用默认的parula色图这两种色图对重金属污染图的观感差异非常大。6.2 用Word COM接口把图插进文档Windows环境下可以用actxserver操作WordWord actxserver(Word.Application); Word.Visible true; doc Word.Documents.Add; for i 1:numel(figFiles) Word.Selection.EndKey(6); % wdStory 移到文档末尾 doc.Content.InsertAfter(sprintf(\n图%d重金属污染分布\n, i)); doc.InlineShapes.AddPicture(fullfile(pwd, figFiles{i})); end doc.SaveAs2(fullfile(pwd, output_report.docx)); doc.Close; Word.Quit;如果用doc.Content.InsertAfter直接插入第二次插入会覆盖文档内容不要这么写必须先写Word.Selection.EndKey(6)把光标移动到文档末尾再插入图片和标题。这里EndKey(6)是Word中的wdStory单位代表整篇文档的结尾。最后一个非常实用的技巧当你在循环里创建多个figure时最好用set(gcf,NumberTitle,off)同时把FigureName设成元素的名称否则figure序号会出现在标题栏影响截图。更重要的是在导出前统一colorbar范围避免8种元素用8个不同的色标干扰看图ax findobj(gcf, Type, axes, Colorbar, on); clim([cmin, cmax]);这个clim在R2022a以后替代了旧的caxis直接写clim更兼容新版Matlab。以上代码片段都可以在真实数据上直接替换使用跑通之后你就能得到一个自动生成的高分辨率图件集连同源代码一起交付。本文还有配套的精品资源点击获取