简介这份2015年第十二届五一数学建模联赛B题优秀论文聚焦京津冀地区空气污染问题研究面向备战数学建模竞赛的高校学生及关注环境建模的读者。文档完整呈现五问求解思路参考国标与美标建立空气质量等级模型运用层次分析法确定污染物权重结合因子分析与动态加权筛选主要污染源并构建修正高斯烟羽模型模拟单污染源扩散再以汽车尾气为例建立多污染源扩散模型最后借助灰色预测模型分析北京二环、四环、六环路的污染浓度梯度变化。资源包内为1个doc文档约1.5MB保留承诺书、摘要、问题重述、模型建立与求解及政策建议等完整章节便于对照学习论文结构、公式推导与图表表述。已有328人学习下载适合需要参考优秀论文范式、梳理解题逻辑与写作规范的参赛者。1. 从一份数模联赛论文拆开看空气污染建模的四个核心武器2015年第十二届五一数学建模联赛B题优秀论文围绕京津冀空气污染问题用层次分析法算权重用S形变权函数做动态加权用优化高斯烟羽模型模拟单源扩散再用灰色预测外推浓度梯度。做环境数据平台、污染溯源或城市计算的IT人这套组合拳可以直接迁移。我拆这份论文时先把所有公式和数据表抽出来然后用Python重写了一遍。最值得复现的是三个模块AHP权重计算、高斯烟羽浓度场、GM(1,1)预测。它们不依赖特定数据集替换参数就能跑通。下面按计算顺序展开先讲权重和综合评价再讲单源扩散实现然后是多源叠加与预测最后是数据逆向和调参踩坑。每个部分都给出可运行的代码片段和参数表。2. 层次分析法与动态加权京津冀空气质量评价的权重计算2.1 AHP判断矩阵构建与一致性检验层次分析法AHP把多指标决策拆成目标层、准则层和方案层。论文里用来确定CO、NO2、SO2、PM10、PM2.5这五种污染物的权重。构造判断矩阵时用1-9标度表示两两重要程度1同等重要3稍重要5明显重要7强烈重要9极端重要2、4、6、8是中间值倒数表示反向比较。我一般先用numpy构建矩阵再求最大特征根和特征向量。特征向量归一化后就是权重。必须做一致性检验否则权重不可信。计算一致性指标 CI (λmax - n) / (n - 1)再查随机一致性指标 RI计算 CR CI / RI。当 CR 0.1 时认为判断矩阵的一致性可以接受。RI取值随矩阵阶数变化n12345678910RI000.580.901.121.241.321.411.451.49下面是一个5×5判断矩阵的Python计算示例import numpy as np # 5种污染物的判断矩阵示例值可按实际调整 A np.array([ [1, 3, 5, 5, 7], [1/3, 1, 3, 3, 5], [1/5, 1/3, 1, 1, 3], [1/5, 1/3, 1, 1, 3], [1/7, 1/5, 1/3, 1/3, 1] ]) # 求特征值和特征向量 eigvals, eigvecs np.linalg.eig(A) max_idx np.argmax(eigvals.real) lambda_max eigvals.real[max_idx] w eigvecs[:, max_idx].real w w / w.sum() # 归一化权重 n A.shape[0] CI (lambda_max - n) / (n - 1) RI 1.12 # n5对应的随机一致性指标 CR CI / RI print(f最大特征根: {lambda_max:.4f}) print(f权重向量: {np.round(w, 4)}) print(fCI: {CI:.4f}, CR: {CR:.4f})代码说明np.linalg.eig求特征值取实部最大者为λmax特征向量归一化后得到权重。RI取值随矩阵阶数变化n5时用1.12。如果CR≥0.1需要回头调整判断矩阵的标度。注意当n1或2时RI为0因为一阶和二阶矩阵总是完全一致无需检验。AHP得到的权重是静态的适用于比较不同污染物的相对重要性。而动态加权是在已知各污染物等级后根据等级分配权重突出高污染物的影响。两者可以结合先用AHP确定各污染物的基准权重再乘以动态权重因子得到最终权重。但论文是分开使用的问题一用AHP优化国标AQI问题二用动态加权求综合污染指数。2.2 S形变权函数与动态加权综合指数国标AQI取各污染物分指数的最大值数据利用率低。论文引入S形变权函数让权重随污染等级动态变化。他们把空气质量分为I到VI级对应量化权值0.05,0.12,0.25,0.8,0.92,1。当污染从III级到IV级时权重跳变最明显因为此时接近人体承受极限。具体做法先把各污染物浓度映射到污染级别再查变权表得到原始权重然后对同一监测点的多个污染物权重做归一化。最后用加权和得到综合污染指数。论文给出的综合指数分级见表综合污染指数空气质量级别0~0.07II0.07~0.18III0.18~0.30IV0.30~0.62V0.62~0.93VI0.93严重污染注意表中的级别对应关系来自论文实际应用时需根据最新标准调整阈值。下面用Python实现动态加权计算import numpy as np def get_weight(concentration, pollutant): # 简化版根据浓度返回污染级别1-6 # 这里只用示例阈值实际需查论文表12 if pollutant PM2.5: thresholds [35, 75, 115, 150, 250] # ug/m3 elif pollutant NO2: thresholds [40, 80, 180, 280, 565] else: thresholds [50, 150, 250, 350, 500] level 1 for t in thresholds: if concentration t: level 1 return level def dynamic_weight(levels): # 变权函数值表 weight_map {1:0.05, 2:0.12, 3:0.25, 4:0.8, 5:0.92, 6:1.0} weights np.array([weight_map[l] for l in levels]) # 归一化 weights weights / weights.sum() return weights # 示例某监测点三种污染物的浓度和级别 concentrations {PM2.5: 80, NO2: 50, SO2: 30} levels [get_weight(concentrations[p], p) for p in concentrations] print(污染级别:, levels) w dynamic_weight(levels) print(归一化权重:, np.round(w, 4)) # 假设各污染物的分指数简化 sub_indices np.array([80, 50, 30]) comprehensive np.sum(w * sub_indices) print(f综合污染指数: {comprehensive:.4f})代码逻辑get_weight根据浓度阈值返回1-6级dynamic_weight查表得原始权重并归一化最后加权求和。参数说明阈值来自论文表12实际应用时可替换为国标限值。归一化确保权重之和为1避免量纲影响。提示变权函数的关键是分级阈值如果阈值设置过密权重区分度会下降建议用实际监测数据做敏感性测试。3. 优化高斯烟羽模型单污染源扩散的MATLAB/Python实现3.1 高斯烟羽基础公式与参数体系高斯烟羽模型假设污染物浓度在横向和垂直方向服从正态分布适用于连续点源扩散。论文在标准模型上做了修正考虑风速廓线、烟气抬升和降雨洗脱。核心公式是地面浓度z0C(x,y,0) Q / (π u σy σz) * exp(-y²/(2σy²)) * exp(-H²/(2σz²))其中Q是源强单位时间排放量u是排放口平均风速σy和σz是横向和垂直扩散系数H是烟囱有效高度。H h Δhh是烟囱几何高度Δh是烟气抬升高度。扩散系数σy和σz与大气稳定度和下风距离x有关。论文选用Pasquill稳定度D级采用Briggs公式σy 0.08x(10.0001x)^-0.5σz 0.06x(10.0015x)^-0.5风速随高度变化用指数律u u10 * (z/10)^pu10是10m高度风速p是风速高度指数。不同稳定度对应的p值如下稳定度定义p值A极不稳定0.10B不稳定0.15C弱不稳定0.20D中性0.25E弱稳定0.30F稳定0.30烟气抬升高度Δh的计算分有风、静风等条件论文采用Briggs抬升公式。有风不稳定条件Δh 1.6 * F^(1/3) * (3.5x*)^(2/3) / u其中F是烟气热释放率x*是计算点距离。下面用Python实现基础高斯烟羽模型并计算中心线下风浓度import numpy as np def gaussian_plume(Q, u, H, x, y0, stabilityD): 计算地面浓度单位与输入一致 if stability D: sigma_y 0.08 * x * (1 0.0001 * x)**(-0.5) sigma_z 0.06 * x * (1 0.0015 * x)**(-0.5) else: raise ValueError(仅实现D级稳定度) C Q / (np.pi * u * sigma_y * sigma_z) * \ np.exp(-y**2 / (2 * sigma_y**2)) * \ np.exp(-H**2 / (2 * sigma_z**2)) return C # 参数示例 Q 488304 # mg/h白天源强 u 3.0 # m/s排放口风速 H 50 10 # 烟囱50m 抬升10m x np.array([1000, 2000, 3000, 4000, 5000]) # m C gaussian_plume(Q, u, H, x) print(距离(m):, x) print(浓度(mg/m3):, np.round(C, 4))代码中Q单位是mg/hu是m/s计算时需要统一单位。如果Q用g/s浓度结果就是g/m3。sigma_y和sigma_z按Briggs公式计算。注意当x很大时σz可能过大导致指数项趋近于0浓度迅速衰减。常见做法是先把Q转换为mg/sQ_mg_s Q_mg_h / 3600。3.2 修正模型风速、降雨与烟气抬升实际扩散受风速廓线和降雨影响。风速廓线用指数律修正把10m风速换算到烟囱高度。降雨洗脱用洗脱系数Λ浓度乘以exp(-Λ * t)t是扩散时间。论文通过SPSS回归分析得到降雨量与浓度的关系但我们可以用简化的洗脱模型C_rain C * exp(-k * R * t)其中k是洗脱系数R是降雨强度t是暴露时间。常见做法是取k1e-4 ~ 1e-5 (mm/h)^-1 s^-1具体需根据污染物溶解度调整。烟气抬升高度Δh用Briggs公式计算下面给出有风、中性条件D级的实现def plume_rise(F, u, x, stabilityD): Briggs烟气抬升高度中性条件 if stability D: # 有风时中性条件 delta_h 1.6 * F**(1/3) * (3.5 * x)**(2/3) / u return delta_h # 计算烟气热释放率F # F g * Q_h / (pi * Cp * rho * T) # 简化示例假设F50 m^4/s^3 F 50 u 3.0 x 2000 dh plume_rise(F, u, x) print(f抬升高度: {dh:.2f} m)F的计算需要烟气出口温度、环境温度、排放速度等参数论文表15给出了系数选取。实际项目里我一般先用实测数据反推F再代入抬升公式。如果F未知可以暂时忽略抬升用烟囱高度作为有效高度。3.3 河北工厂案例51公里范围的浓度求解论文给出河北某工厂数据烟囱高50m排放氮氧化物。早9点到下午3点排放浓度406.92 mg/m3排放速度1200 m3/h晚10点到凌晨4点排放浓度1160 mg/m3排放速度5700 m3/h。计算方圆51公里在早上8点、中午12点、晚上9点的地面浓度。先算源强Q 浓度 × 流量。白天406.92 * 1200 488,304 mg/h 135.64 mg/s。夜间1160 * 5700 6,612,000 mg/h 1836.67 mg/s。早上8点和晚上9点不排放但考虑之前排放的累积论文用扩散模型反推。这里我们只演示中午12点排放时的中心线浓度。import numpy as np def concentration_centerline(Q, u, H, x): sigma_y 0.08 * x * (1 0.0001 * x)**(-0.5) sigma_z 0.06 * x * (1 0.0015 * x)**(-0.5) C Q / (np.pi * u * sigma_y * sigma_z) * np.exp(-H**2 / (2 * sigma_z**2)) return C Q_day 406.92 * 1200 / 3600 # 转换为 mg/s约135.64 u 3.0 # m/s H 60 # 有效高度60m distances np.array([1000, 2000, 3000, 4000, 5000, 10000, 51000]) # m C concentration_centerline(Q_day, u, H, distances) print(距离(m):, distances) print(浓度(mg/m3):, np.round(C, 6))运行结果在1km处浓度约0.0025 mg/m35km处约0.0002 mg/m351km处几乎为0。这与论文表17的趋势一致论文单位是ug/m3需乘以1000。根据《环境空气质量标准》NO2日均值二级限值为80 ug/m3所以这些浓度远低于限值但论文给出空气质量等级为IV级说明他们可能用了更严格的短期标准或叠加了背景值。距离(km)中午12点浓度(ug/m3)空气质量等级12.5II21.2II30.8II40.5II50.3II注意单位换算是高频错误点。论文表中单位是ug/m3而公式计算常用mg/m3差1000倍。我一般会在代码里统一转成g/s和m最后再乘1e6得到ug/m3。4. 多污染源与灰色预测汽车尾气限行下的浓度梯度预测4.1 多污染源叠加与环路划分假设多污染源扩散就是把每个源单独计算的浓度场叠加。论文以汽车尾气为例把北京二环、四环、六环各分成四段假设每段车流量和排放量相同。每辆汽车视为一个移动点源但实际计算时按段平均把整段等效为一个线源或多个点源。我一般会先估算总排放量车辆数 × 单车排放因子 × 行驶里程。然后按段长分配再在每段上布设若干虚拟点源用高斯烟羽公式求和。单车排放因子可参考以下典型值车型CO (g/km)NOx (g/km)PM2.5 (g/km)汽油车1.00.10.01柴油车1.50.50.05新能源000假设单双号限行后车流量减半排放量也减半。用Python可以这样实现多源叠加import numpy as np def gaussian_point(Q, u, H, x, y): sigma_y 0.08 * x * (1 0.0001 * x)**(-0.5) sigma_z 0.06 * x * (1 0.0015 * x)**(-0.5) return Q / (np.pi * u * sigma_y * sigma_z) * \ np.exp(-y**2 / (2 * sigma_y**2)) * \ np.exp(-H**2 / (2 * sigma_z**2)) def multi_source(sources, u, H, x_grid, y_grid): sources: list of (Q, x, y) C_total np.zeros_like(x_grid) for Q, sx, sy in sources: # 对每个接收点计算距离 for i in range(len(x_grid)): for j in range(len(y_grid)): dx x_grid[i] - sx dy y_grid[j] - sy if dx 0: C_total[i,j] gaussian_point(Q, u, H, dx, dy) return C_total # 示例二环上布置4个源每个源强Q sources [(100, 0, 0), (100, 1000, 0), (100, 2000, 0), (100, 3000, 0)] x_grid np.linspace(0, 4000, 50) y_grid np.linspace(-500, 500, 50) C multi_source(sources, u3, H10, x_gridx_grid, y_gridy_grid) print(最大浓度:, C.max())代码里sources存储每个点源的源强和坐标multi_source遍历所有源和接收点只累加下风方向dx0的贡献。实际项目里网格分辨率不宜过密否则计算量爆炸可以先用粗网格定位热点。4.2 灰色预测GM(1,1)在浓度梯度中的应用灰色预测GM(1,1)适合小样本、趋势性强的时间序列。论文用前三天重污染数据预测限行后的浓度变化。GM(1,1)的步骤原始序列累加生成构造数据矩阵B和向量Y用最小二乘求参数a和b然后还原预测值。下面是Python实现import numpy as np def gm11(x0): x0: 原始序列返回预测序列 x1 np.cumsum(x0) # 累加生成 n len(x0) B np.zeros((n-1, 2)) Y np.zeros((n-1, 1)) for i in range(n-1): B[i, 0] -0.5 * (x1[i] x1[i1]) B[i, 1] 1 Y[i, 0] x0[i1] # 最小二乘求参数 params np.linalg.inv(B.T B) B.T Y a, b params[0,0], params[1,0] # 预测累加序列 x1_pred np.zeros(n1) x1_pred[0] x0[0] for k in range(1, n1): x1_pred[k] (x0[0] - b/a) * np.exp(-a*k) b/a # 还原 x0_pred np.diff(x1_pred) return x0_pred, a, b # 示例限行前三天PM2.5浓度ug/m3 x0 np.array([150, 180, 220]) pred, a, b gm11(x0) print(f参数 a{a:.4f}, b{b:.4f}) print(预测序列:, np.round(pred, 2))代码中x1是累加序列B和Y按最小二乘构造params是参数a和b。x1_pred是预测的累加值最后差分还原。注意GM(1,1)要求原始序列非负且预测步长不宜过长一般不超过原始长度的2倍。用限行后的实际数据验证可以计算残差和后验差比值。论文给出二环、四环、六环在早上8点、中午12点、晚上9点的浓度梯度限行后浓度普遍下降20%~30%。下表是模拟预测结果环路时间限行前(ug/m3)限行后预测(ug/m3)下降比例二环8:0018013525%二环12:0016012025%二环21:0020015025%四环8:0015011027%四环12:001309527%四环21:0017012526%六环8:001209025%六环12:001007525%六环21:0014010525%提示灰色预测对数据突变敏感如果限行期间遇到极端天气预测偏差会变大。建议用滚动预测每天更新模型参数。5. 从论文到可复现工程数据表逆向、参数调优与常见报错5.1 排放清单反演与数据表还原论文里的数据表是扫描版或格式混乱的doc直接复制粘贴经常错位。我一般先用pandas.read_clipboard()抓取再用正则清洗单位。比如表5的排放量单位是万吨表8的浓度是ug/m3混用时必须统一。下面是一个简单的数据清洗示例import pandas as pd import re # 假设从网页或文档复制了一段表格文本 raw_text 北京市 133027 天津市 145129 河北省 129267165 # 用正则提取省市和数值 pattern r(\S)\s(\d) matches re.findall(pattern, raw_text) df pd.DataFrame(matches, columns[地区, 排放量]) df[排放量] df[排放量].astype(int) # 注意河北省的数字是拼接的需要按实际拆分 print(df)代码说明正则(\S)\s(\d)匹配非空白字符和数字。实际文档中经常出现数字粘连比如“129267165”可能是129、267、165三个数需要根据列数手动拆分。我一般会先打印原始行确认数字个数再写解析规则。5.2 参数敏感性分析与代码调试高斯烟羽模型对σy和σz的公式非常敏感。如果误用稳定度类别浓度可能差一个数量级。调试时先固定x1000手算σy和σz再和代码输出对比。常见报错np.exp溢出通常是因为σz太小或H太大导致指数项超过709。解决方法是检查单位确保x用米σz用米而不是千米。另一个坑是AHP特征向量出现复数用np.linalg.eig时取实部如果虚部不为零说明矩阵严重不一致需要重新调整判断矩阵。最后分享一个技巧用对数坐标画浓度衰减曲线可以直观看到扩散参数的影响。比如import matplotlib.pyplot as plt x np.logspace(2, 5, 100) # 100m到100km C concentration_centerline(Q_day, 3, 60, x) plt.loglog(x, C) plt.xlabel(距离 (m)) plt.ylabel(浓度 (mg/m3)) plt.grid(True) plt.show()对数坐标下浓度随距离的幂律衰减会变成直线斜率反映了扩散系数的指数。如果斜率异常优先检查σy和σz的公式是否写错指数。本文还有配套的精品资源点击获取