简介一份面向环境科学与气候政策研究的MATLAB建模文档系统呈现温室气体排放模型的构建与应用路径。文档首先梳理工业、交通、农业、能源生产等主要温室气体排放源的识别方法以及二氧化碳、甲烷、氧化亚氮等数据获取渠道随后基于能量守恒原理建立分层大气能量平衡方程借助矩阵运算、差分方法与数值积分求解辐射传输过程通过迭代直到能量平衡模拟温度场随温室气体浓度变化的趋势并与实际观测数据进行校验。在减排策略层面模型可引入碳税税率、清洁能源使用比例等参数定量比较不同政策情景下温室气体浓度的变化为科学决策提供依据。资源包含一个docx文档压缩包整体约34KB共1个文件内容精炼但纵向覆盖建模全流程。适合具备一定MATLAB操作能力的环境科学家、气候研究人员、政策制定者与开发者研读也可作为高校相关课程的案例参考。目前已有61人学习阅读后能形成从排放数据整理、数值仿真到减排策略评估的完整分析框架。1. 温室气体排放建模为什么第一步是算清清单而不是直接调优化温室气体排放建模与减排策略分析第一步往往不是急着写优化代码而是把排放清单先算明白。你拿到的原始资料大概率是几十张 Excel分燃料的消费量、外购电量、产品产量、运输里程配上来源各异的排放因子清单一旦对不上账后面所有情景推演和策略寻优都在错误基础上空转。MATLAB 在这个场景里合适的理由很直接核算、可视化、情景推演、优化求解可以在同一个工作空间里闭环不用在 Excel、Python 和报告之间来回搬运数据。这套方案不需要特殊硬件也不需要额外数据源只要核算表齐全就能落地适合做能源环境数据分析的工程师、参加数学建模的同学以及要给业务方解释减排策略是怎么算出来的分析师。下面按核算、推演、寻优、验证四个环节展开每一段都给出能直接抄走改用的代码。2. 用排放因子法在MATLAB里搭建温室气体排放核算模型排放因子法IPCC 清单指南里最常用的一层核算方法的核心公式只有一行排放量等于活动数据乘以排放因子。活动数据回答做了多少比如烧了多少吨煤、用了多少兆瓦时电、生产了多少吨水泥排放因子回答每做一单位排放多少来自实测、行业默认值或国际缺省值。公式朴素但真实的建模工作量全在数据口径上核算边界划到哪一层活动数据和排放因子用什么字段对齐CH4、N2O 怎么折算成 CO2 当量。口径没定模型代码写得再漂亮都是白搭。2.1 核算边界与数据表先定口径再写代码常见的做法是把核算对象分成四类化石燃料燃烧、外购电力和热力、工业生产过程、废弃物处理。每一类对应一张活动数据表和一张排放因子表两张表通过一个唯一编号关联。这样做的原因是活动数据的统计口径差异很大燃料台账是按实物量记的电力是按兆瓦时记的水泥产量是按吨记的只有先统一成表 主键的结构后续才能用 MATLAB 的 table 直接做集合运算。排放源类别典型活动数据排放因子单位常用数据来源化石燃料燃烧分燃料消费量吨或万Nm³tCO2/t 燃料tCO2/万Nm³企业能源台账、区域能源平衡表外购电力/热力年度用电量MWhtCO2/MWh电网排放因子、供热企业报告工业生产过程产品产量ttCO2/t 产品行业默认因子水泥、钢铁等废弃物处理填埋量、焚烧量ttCO2/t 废物IPCC 默认因子然后是 GWP 折算。CH4 和 N2O 需要按 100 年口径折算成 CO2 当量CO2e按 IPCC 第六次评估报告的常用取值CH4 取 28、N2O 取 273。这个选择直接影响总量写报告时必须注明用的是哪套口径否则第二年跟别人对数据时对不上。我的习惯是在因子表里放一列 GWP和排放因子分开存核算时再乘这样改口径只需要改一列。2.2 用 table 组织活动数据与排放因子避免索引错位的三个习惯第一表里放一个唯一主键我一般用 src_id不要靠表内顺序对齐。Excel 表经过几手整理后行序经常变化按顺序相乘是最隐蔽的错误来源。第二readtable 读 Excel 时加 VariableNamingRule, preserve否则中文列名或者带空格的列名会被自动改成 Var1、Var2后面所有引用全乱。第三单位在进入计算前统一折算成 tCO2e不要把万吨和吨混在一个公式里宁可多写一行注释也不要在心里记换算。% 读取活动数据和排放因子表保留中文列名 activity readtable(activity_data.xlsx, Sheet, activity, ... VariableNamingRule, preserve); ef readtable(ef.xlsx, Sheet, ef, ... VariableNamingRule, preserve); % 按 src_id 对齐两张表防止行序错位 aggr outerjoin(activity, ef(:, {src_id, EF, GWP}), ... Keys, src_id, MergeKeys, true); % 核心核算活动量 × 排放因子 × GWP单位统一到 tCO2e/yr aggr.E_t aggr.Activity .* aggr.EF .* aggr.GWP; % 按行业汇总并输出 result groupsummary(aggr, sector, sum, E_t); writetable(result, emission_result.csv);这段代码的关键在 outerjoin它按 Keys 指定的 src_id 做外连接MergeKeys, true 表示两个表里同名的 src_id 合并成一列而不是各自保留。活动数据表里除了 src_id 和 Activity 之外还有年份、行业等列它们会被原样带进 aggr。如果两张表里有其他同名列MATLAB 会自动加 _activity、_ef 这样的后缀引用前先用 whos 确认列名。注意outerjoin 对同名列自动加后缀下次引用时用 aggr.EF 之前先 whos aggr看清列名再写避免变量名拼写错误在后期才暴露。核算那行的.*是关键操作符。MATLAB 对数组做逐元素运算必须用点乘写成aggr.Activity * aggr.EF会触发矩阵乘法报错或者算出完全错误的结果。对几千行的季度数据来说矢量化比 for 循环快一个数量级而且代码短一半。groupsummary 的第一个参数是表第二个是分组列 sector第三和第四个参数分别是聚合方法与聚合列返回的表里自动生成 sum_E_t 列名后续引用要按这个新列名来。2.3 矢量化核算代码一行算完一年的排放往上面这段里再加一个维度就是年度核算。活动数据表里通常有 year 列先按年份和行业排好序再用 reshape 把长表转成宽表就能直接画堆叠图。排序这步不能省reshape 是按列填充的如果数据没排序矩阵的行列对应关系就是乱的画出来的图看似正常实际上年份和行业全错位了。% 先按年份、行业排序保证 reshape 顺序一致 result sortrows(aggr, {year, sector}); nYear 10; nSector 4; % 重排成 年份×行业 矩阵每行一年每列一个行业 Emat reshape(result.E_t, nYear, nSector); % 堆叠柱状图柱段累加适合看结构变化 figure; bar(Emat, stacked); legend(sectorNames, Location, northwest); xlabel(年份); ylabel(排放量 (万tCO2e));reshape 按列填充所以排序时 sector 必须在 year 之后作为二级排序键保证每个年份内部的行业顺序固定。sectorNames 是行业的元胞数组顺序必须和 Emat 的列顺序一致这个对应关系交稿前值得人工核对一次。如果数据是逐月的把 nYear 换成 nMonth、把年份标签换成月份标签即可代码结构不用动。2.4 排放结构可视化从核算表到能直接进报告的图堆叠柱状图适合看结构变迁看具体的行业强弱排序则用 heatmap。MATLAB 的 heatmap 直接吃矩阵把 年份×行业 的 Emat 传进去就能出结果横轴显示年份、纵轴显示行业颜色深浅代表排放量高低。报告中常需要哪一年哪个行业贡献最大这种结论heatmap 一眼就能读出来。如果还要对比不同情景把三张 Emat 放进同一张 figure 用 tiledlayout 并排展示比反复翻 Excel 高效得多。% 用热力图检查是否存在异常突变的年份 figure; heatmap(yearList, sectorNames, Emat); colormap(parula); colorbar;热力图在这里还有一个用途排错。某个行业某年数据如果异常偏大或偏小heatmap 上会出现明显色块比看数字表更容易发现问题。常见的原因包括活动数据录入单位错误、排放因子这一年被替换、以及某年缺数据被补了 0。发现异常后回头查原始表比对着公式找错快得多。3. 减排策略情景推演在MATLAB里参数化所有如果减排策略本身并不直接减它是通过改变活动水平和排放因子间接改变排放的。提高能效单位产值能耗下降活动数据变小电源结构清洁化单位电量的排放因子变小电气化率上升终端化石燃料被电力替代活动结构发生迁移。把这些策略翻译成参数模型就变成了一台可以反复按下如果键的推演器。关键在参数要少而全参数太多组合空间爆炸分析无从下手参数太少又说不清不同策略之间的差别。3.1 把减排措施翻译成可调参数的情景参数表我一般建一个 scenario 参数表而不是把参数散落在脚本里。散落变量的问题是改一处漏一处表结构让参数、取值和说明待在一起后面按行循环和输出对比都干净。基准情景取最近三年的实际均值政策情景取行业规划或可研报告里的目标值每个参数都要能说清来源。参数物理含义基准情景减排情景数据依据alpha能源强度年下降率%0.51.2行业节能目标、技改报告renew可再生能源发电占比%2035电源发展规划elecShare电力排放占总排放比例%3535由第 2 章核算结果计算beta需求侧增长不确定性%22产值预测区间elecShare 不是策略变量它是从核算模型里算出来的结构参数固定下来减小分析维度。alpha 和 renew 才是策略旋钮后面做交叉推演时只动这两个。beta 放进蒙特卡洛而不是情景表因为它描述的是不确定性而不是策略选择二者混在一起会让结果解读变得困难。3.2 用 meshgrid 展开多策略组合9 个情景一次算完两个参数各取三个水平全部组合是 9 个情景。meshgrid 是展开这种组合最直接的工具生成两个矩阵经过 I(:) 展平后拼进一张表每行一个情景。三个以上参数做交叉时组合数长得很快所以一般只对最关心的 2 到 3 个参数做交叉其余参数固定为中间值。% 两个策略参数各取 3 个水平共 9 个组合情景 alphaList [0.5, 0.8, 1.2]; % 能源强度下降率 renewList [20, 28, 35]; % 可再生能源占比 [A, R] meshgrid(alphaList, renewList); scen table(A(:), R(:), VariableNames, {alpha, renew}); % 基准年排放(万tCO2e)由第 2 章核算结果直接读入 E0 100; elecShare 0.35; Emat zeros(height(scen), 10); % 逐年推演强度下降是复利电力清洁化按占比折算 for k 1:height(scen) for t 0:9 Emat(k, t1) E0 * (1 - scen.alpha(k)/100)^t ... * (1 - scen.renew(k)/100 * elecShare); end end推演公式里(1 - alpha/100)^t是能效改进的复利效果t 从 0 到 9 表示推演 10 年。(1 - renew/100 * elecShare)表达的是可再生能源占比每提高 1 个百分点大约有 elecShare 比例的排放跟着下降。这里把电力部门和非电力部门合并成一行简化了模型代价是无法单独观察电力行业的路径需要细化时把电力排放单独拆一列而非电力部分完全不受 renew 影响。循环里 Emat 的行对应情景、列对应年份画图时直接plot(Emat)就能得到 9 条曲线组成的曲线族。3.3 蒙特卡洛抽样用拉丁超立方给减排路径画置信带点估计只能给出一条线业务方下一步一定会问如果参数没那么乐观怎么办。把 alpha、renew、beta 当成随机变量抽样跑几百次用分位数区间回答这类问题。抽样方法用 lhsdesign 而不是 rand原因是拉丁超立方把参数空间分层后再取样300 到 500 组样本就能稳定覆盖而普通随机抽样需要更多样本才能达到同样的覆盖率。n 500; % 抽样次数300~500 足够稳定 params lhsdesign(n, 3); % 3 列 [0,1] 均匀分层样本 % 把 [0,1] 映射到参数实际区间 alpha 0.5 params(:,1) * 0.9; % 能效提升率 0.5% ~ 1.4% renew 20 params(:,2) * 15; % 可再生占比 20% ~ 35% beta 0.02 params(:,3) * 0.03;% 需求侧年波动 2% ~ 5% E_sim zeros(n, 10); E0 100; elecShare 0.35; for i 1:n for t 0:9 E_sim(i, t1) E0 * (1 - alpha(i)/100)^t ... * (1 - renew(i)/100 * elecShare) * (1 beta(i))^t; end end % 逐年取 5%、50%、95% 分位数 pct prctile(E_sim, [5, 50, 95], 1);映射公式是min 样本 * (max - min)把 [0,1] 的样本平移到真实区间。prctile 的第三个参数 1 表示沿第一维样本维逐列计算得到 3×10 的结果第一行是最悲观线第二行是中位数第三行是最乐观线。画图时用 fill 填充中位数到上下界的区间带比画几十条散线清楚。另外一个细节beta 带 t 次方意味着需求侧的不确定性会随时间放大这正是蒙特卡洛要暴露出来的信息——长期预测的置信带必然比短期宽宽到什么程度是决策者该看见的。4. 组合BP神经网络与优化工具箱寻优减排策略4.1 为什么核算模型之外还要 BP 神经网络拟合曲线核算模型是白箱每一步都有物理含义但它依赖活动数据和排放因子齐全。现实中常见的情况是只有产值、人口、气温这类间接指标排放和它们的关系非单调、说不清机理这时数据驱动的黑箱拟合就有用。BP 神经网络在 MATLAB 里用 feedforwardnet 就能搭典型用途有两个一是做短期排放预测输入经济指标和气象变量输出年度排放二是做排放因子的反标定用实测序列反推等效因子。要提醒的是 BP 替代不了核算它只做拟合和预测策略比选时也有人用层次分析法AHP打分AHP 适合定权重的决策场合想讲客观最优仍然要把目标函数和约束写清楚交给优化求解更可复现。4.2 MATLAB 里训练 BP 网络的最小流程与参数设置% 输入经济产值、人口、能源消费强度输出年度排放 X [gdp; pop; intensity]; % 3×N 输入矩阵 Y emissions; % 1×N 输出向量 % 归一化到 [0,1]psX/psY 保存映射参数用于反归一化 [Xn, psX] mapminmax(X, 0, 1); [Yn, psY] mapminmax(Y, 0, 1); % 单隐藏层 10 个神经元trainlm 收敛快但吃内存 net feedforwardnet(10, trainlm); net.trainParam.epochs 1000; net.trainParam.goal 1e-5; net.trainParam.min_grad 1e-7; % 时间序列必须按顺序切分最后 20% 做测试 net.divideFcn divideblock; net.divideParam.trainRatio 0.8; net.divideParam.valRatio 0; net.divideParam.testRatio 0.2; % 训练并做全样本回代 [net, tr] train(net, Xn, Yn); Y_pred mapminmax(reverse, net(Xn), psY); test_idx tr.testInd; mape mean(abs((Y_pred(test_idx) - Y(test_idx)) ./ Y(test_idx))) * 100;输入矩阵约定每列一个样本、每行一个特征输出是行向量。mapminmax 把数据压到 [0,1] 避免神经元饱和训练完必须用 ps 结构做反归一化否则预测值和真实值根本不在一个量级。feedforwardnet(10, trainlm) 表示单隐藏层 10 个神经元trainlm 是 Levenberg-Marquardt中等样本收敛快数据噪声大时换成 trainbr 更稳。divideFcn 默认是随机切分时间序列必须改 divideblock把最后两年留作测试否则测试集混在样本中间等于提前偷看了未来信息。提示样本不到 100 个时隐藏层神经元从 5 个起步别一上来就 20。神经元越多越容易在训练集上拟合得漂亮测试集上一塌糊涂。tr 结构里存了训练和测试索引tr.best_epoch 记录最佳迭代轮次作图时可以画出训练误差和测试误差的下降曲线看两者是否在某一轮开始分叉。如果测试集 MAPE 在 5% 以下模型基本可以对外讲结论超过 10% 就要回去检查是输入变量漏了关键驱动因素还是神经元数量不对而不是急着堆数据。4.3 将减排策略寻优写成带约束的优化问题减排策略寻优的标准写法是决策变量 x_i 表示第 i 项措施的实施水平连续取值 0 到 1。目标是总成本最小约束有两个年度排放不超过上限总成本不超过预算。减排量要和核算模型的基准 E0 对应用年口径如果是跨年规划把逐年约束展开每一年的排放都写成 x 的线性函数。这类问题在 MATLAB 优化工具箱里用 optimproblem 描述最省事目标函数和约束写成数值表达式求解器会自动匹配算法。4.4 用 MATLAB 优化工具箱的 linprog 求解并检查解的合理性% 5 项减排措施实施水平 0~1 x optimvar(x, 5, LowerBound, 0, UpperBound, 1); unitCost [12; 8; 25; 40; 6]; % 单位减排成本万元/万tCO2e reduction [30; 18; 55; 90; 10]; % 各措施年减排量万tCO2e E0 100; cap 60; budget 800; % 基准排放、排放上限、预算上限 prob optimproblem(ObjectiveSense, minimize); prob.Objective unitCost * x; % 总成本最小 prob.Constraints.emissionCap E0 - reduction * x cap; prob.Constraints.budget unitCost * x budget; [sol, fval] solve(prob, Options, ... optimoptions(linprog, Display, iter));optimvar 声明决策变量并给出边界0 表示不实施1 表示满额实施。排放约束写成E0 - 减排量 cap的原因是引入措施后排放等于基准减去减排量。solve 检测到目标与约束都是线性的会自动调用 linprog如果目标或约束里出现非线性需要换成 fmincon它要求提供初值且对初值敏感建议跑多个随机起点。连续解 x 等于 0.6 表示实施到六成实际项目可以按此排期如果措施必须整体上马把 optimvar 加 Type, integer求解器会切到混合整数规划。Display, iter 打印每轮迭代信息用来判断求解过程是正常收敛还是在约束边界上反复摆动。问题特征推荐求解器需要注意的地方线性目标 线性约束linprogsolve 自动选择稳定但只能处理线性问题非线性目标或约束fmincon对初值敏感需要多起点重试混合整数/大规模组合ga 或 intlinprog规模大时耗时明显先设小种群验证求解结果先做物理性检查各项措施的水平是否都在 0 到 1 之间总成本是否恰好贴着预算或排放上限。如果某个约束明明没到边界说明这个约束不起作用可以删掉缩小问题规模如果解里某项措施是 0但它旁边有个成本更高、减排更少的选择大概率是约束写漏了回头排查约束矩阵。5. 模型验证与敏感性分析让温室气体排放模型可交付5.1 用时间序列留出法回代验证模型模型交付前必须回答凭什么信你。核算模型的验证方式是拿历史年份的活动数据回代把算出的排放总量和公开统计口径做对比误差在几个百分点内说明因子选取合理。BP 模型则靠测试集说话时间序列数据不能随机打乱切分必须按顺序留出最后两年代码里 divideFcn 设为 divideblock 就是这个目的。测试集 MAPE 低于 5% 可以对外讲结论10% 以上就是模型结构有问题继续加数据往往没用应该回去检查输入变量是否漏了关键驱动因素。5.2 局部敏感性分析找出最该优先采集的数据局部敏感性分析的做法是逐个参数上下浮动 10%观察输出变化率。归一化敏感度大于 1 说明输出被参数放大这类参数就是数据采集的优先级所在——它的误差对结果影响最大。把第 3 章的推演逻辑封装成一个函数 run_sim敏感分析、情景推演、优化寻优共用同一个模型避免优化用的是新版本、敏感性用的还是旧版本这类交付事故。params0 [0.8, 25, 0.3]; % 基准参数alpha、renew、beta delta 0.1; sens zeros(numel(params0), 1); for i 1:numel(params0) p_up params0; p_up(i) p_up(i) * (1 delta); p_dn params0; p_dn(i) p_dn(i) * (1 - delta); E_up run_sim(p_up); % 返回 10 年后排放 E_dn run_sim(p_dn); % 归一化敏感度输出相对变化 / 输入相对变化 sens(i) ((E_up - E_dn) / mean([E_up, E_dn])) / (2 * delta); end figure; barh(sens); set(gca, YTickLabel, {alpha, renew, beta}); xlabel(归一化敏感度);分母里的 2 是上下浮动造成的总相对变化这样算出的敏感度就是单位相对变化下的输出响应三个参数之间可以互相比较。alpha 和 renew 通常是最敏感的两个因为它们直接进入复利公式beta 的影响随时间放大所以推演年份越长它的敏感度排名越高。注意局部敏感性分析要求输出对参数连续。参数是离散取值的场景比如月份取整、措施只能整年上马不适合用这个方法改用全局敏感性分析逐个组合跑一遍看分布。5.3 交付前必查的四个 MATLAB 数值坑第一单位。活动数据是吨还是万吨排放因子是 tCO2/t 还是 gCO2/kWh混在一起算出来的数字能差好几个数量级。我的做法是在核算脚本开头统一一个注释块把所有单位标注清楚再写公式。第二GWP 口径。100 年和 20 年基准下 CH4 的折算值差了快三倍报告里写明了哪套口径别人才能复现你的数字。第三中文列名。readtable 不加 VariableNamingRule,preserve 时中文列名会被改成 Var1、Var2等脚本跑完才发现引用全错浪费的时间足够把代码重写一遍。第四优化初值。fmincon 换个 x0 结果可能完全不同至少跑 5 个随机起点或者先用 ga 出一组接近最优的解作为初值再交给 fmincon 精修。 交付时把敏感性分析结果和核算口径说明打包给业务方让对方标出哪些数据最没底这一步比模型精度更能建立信任。本文还有配套的精品资源点击获取 MATLAB温室气体排放建模与减排策略优化实战指南
简介一份面向环境科学与气候政策研究的MATLAB建模文档系统呈现温室气体排放模型的构建与应用路径。文档首先梳理工业、交通、农业、能源生产等主要温室气体排放源的识别方法以及二氧化碳、甲烷、氧化亚氮等数据获取渠道随后基于能量守恒原理建立分层大气能量平衡方程借助矩阵运算、差分方法与数值积分求解辐射传输过程通过迭代直到能量平衡模拟温度场随温室气体浓度变化的趋势并与实际观测数据进行校验。在减排策略层面模型可引入碳税税率、清洁能源使用比例等参数定量比较不同政策情景下温室气体浓度的变化为科学决策提供依据。资源包含一个docx文档压缩包整体约34KB共1个文件内容精炼但纵向覆盖建模全流程。适合具备一定MATLAB操作能力的环境科学家、气候研究人员、政策制定者与开发者研读也可作为高校相关课程的案例参考。目前已有61人学习阅读后能形成从排放数据整理、数值仿真到减排策略评估的完整分析框架。1. 温室气体排放建模为什么第一步是算清清单而不是直接调优化温室气体排放建模与减排策略分析第一步往往不是急着写优化代码而是把排放清单先算明白。你拿到的原始资料大概率是几十张 Excel分燃料的消费量、外购电量、产品产量、运输里程配上来源各异的排放因子清单一旦对不上账后面所有情景推演和策略寻优都在错误基础上空转。MATLAB 在这个场景里合适的理由很直接核算、可视化、情景推演、优化求解可以在同一个工作空间里闭环不用在 Excel、Python 和报告之间来回搬运数据。这套方案不需要特殊硬件也不需要额外数据源只要核算表齐全就能落地适合做能源环境数据分析的工程师、参加数学建模的同学以及要给业务方解释减排策略是怎么算出来的分析师。下面按核算、推演、寻优、验证四个环节展开每一段都给出能直接抄走改用的代码。2. 用排放因子法在MATLAB里搭建温室气体排放核算模型排放因子法IPCC 清单指南里最常用的一层核算方法的核心公式只有一行排放量等于活动数据乘以排放因子。活动数据回答做了多少比如烧了多少吨煤、用了多少兆瓦时电、生产了多少吨水泥排放因子回答每做一单位排放多少来自实测、行业默认值或国际缺省值。公式朴素但真实的建模工作量全在数据口径上核算边界划到哪一层活动数据和排放因子用什么字段对齐CH4、N2O 怎么折算成 CO2 当量。口径没定模型代码写得再漂亮都是白搭。2.1 核算边界与数据表先定口径再写代码常见的做法是把核算对象分成四类化石燃料燃烧、外购电力和热力、工业生产过程、废弃物处理。每一类对应一张活动数据表和一张排放因子表两张表通过一个唯一编号关联。这样做的原因是活动数据的统计口径差异很大燃料台账是按实物量记的电力是按兆瓦时记的水泥产量是按吨记的只有先统一成表 主键的结构后续才能用 MATLAB 的 table 直接做集合运算。排放源类别典型活动数据排放因子单位常用数据来源化石燃料燃烧分燃料消费量吨或万Nm³tCO2/t 燃料tCO2/万Nm³企业能源台账、区域能源平衡表外购电力/热力年度用电量MWhtCO2/MWh电网排放因子、供热企业报告工业生产过程产品产量ttCO2/t 产品行业默认因子水泥、钢铁等废弃物处理填埋量、焚烧量ttCO2/t 废物IPCC 默认因子然后是 GWP 折算。CH4 和 N2O 需要按 100 年口径折算成 CO2 当量CO2e按 IPCC 第六次评估报告的常用取值CH4 取 28、N2O 取 273。这个选择直接影响总量写报告时必须注明用的是哪套口径否则第二年跟别人对数据时对不上。我的习惯是在因子表里放一列 GWP和排放因子分开存核算时再乘这样改口径只需要改一列。2.2 用 table 组织活动数据与排放因子避免索引错位的三个习惯第一表里放一个唯一主键我一般用 src_id不要靠表内顺序对齐。Excel 表经过几手整理后行序经常变化按顺序相乘是最隐蔽的错误来源。第二readtable 读 Excel 时加 VariableNamingRule, preserve否则中文列名或者带空格的列名会被自动改成 Var1、Var2后面所有引用全乱。第三单位在进入计算前统一折算成 tCO2e不要把万吨和吨混在一个公式里宁可多写一行注释也不要在心里记换算。% 读取活动数据和排放因子表保留中文列名 activity readtable(activity_data.xlsx, Sheet, activity, ... VariableNamingRule, preserve); ef readtable(ef.xlsx, Sheet, ef, ... VariableNamingRule, preserve); % 按 src_id 对齐两张表防止行序错位 aggr outerjoin(activity, ef(:, {src_id, EF, GWP}), ... Keys, src_id, MergeKeys, true); % 核心核算活动量 × 排放因子 × GWP单位统一到 tCO2e/yr aggr.E_t aggr.Activity .* aggr.EF .* aggr.GWP; % 按行业汇总并输出 result groupsummary(aggr, sector, sum, E_t); writetable(result, emission_result.csv);这段代码的关键在 outerjoin它按 Keys 指定的 src_id 做外连接MergeKeys, true 表示两个表里同名的 src_id 合并成一列而不是各自保留。活动数据表里除了 src_id 和 Activity 之外还有年份、行业等列它们会被原样带进 aggr。如果两张表里有其他同名列MATLAB 会自动加 _activity、_ef 这样的后缀引用前先用 whos 确认列名。注意outerjoin 对同名列自动加后缀下次引用时用 aggr.EF 之前先 whos aggr看清列名再写避免变量名拼写错误在后期才暴露。核算那行的.*是关键操作符。MATLAB 对数组做逐元素运算必须用点乘写成aggr.Activity * aggr.EF会触发矩阵乘法报错或者算出完全错误的结果。对几千行的季度数据来说矢量化比 for 循环快一个数量级而且代码短一半。groupsummary 的第一个参数是表第二个是分组列 sector第三和第四个参数分别是聚合方法与聚合列返回的表里自动生成 sum_E_t 列名后续引用要按这个新列名来。2.3 矢量化核算代码一行算完一年的排放往上面这段里再加一个维度就是年度核算。活动数据表里通常有 year 列先按年份和行业排好序再用 reshape 把长表转成宽表就能直接画堆叠图。排序这步不能省reshape 是按列填充的如果数据没排序矩阵的行列对应关系就是乱的画出来的图看似正常实际上年份和行业全错位了。% 先按年份、行业排序保证 reshape 顺序一致 result sortrows(aggr, {year, sector}); nYear 10; nSector 4; % 重排成 年份×行业 矩阵每行一年每列一个行业 Emat reshape(result.E_t, nYear, nSector); % 堆叠柱状图柱段累加适合看结构变化 figure; bar(Emat, stacked); legend(sectorNames, Location, northwest); xlabel(年份); ylabel(排放量 (万tCO2e));reshape 按列填充所以排序时 sector 必须在 year 之后作为二级排序键保证每个年份内部的行业顺序固定。sectorNames 是行业的元胞数组顺序必须和 Emat 的列顺序一致这个对应关系交稿前值得人工核对一次。如果数据是逐月的把 nYear 换成 nMonth、把年份标签换成月份标签即可代码结构不用动。2.4 排放结构可视化从核算表到能直接进报告的图堆叠柱状图适合看结构变迁看具体的行业强弱排序则用 heatmap。MATLAB 的 heatmap 直接吃矩阵把 年份×行业 的 Emat 传进去就能出结果横轴显示年份、纵轴显示行业颜色深浅代表排放量高低。报告中常需要哪一年哪个行业贡献最大这种结论heatmap 一眼就能读出来。如果还要对比不同情景把三张 Emat 放进同一张 figure 用 tiledlayout 并排展示比反复翻 Excel 高效得多。% 用热力图检查是否存在异常突变的年份 figure; heatmap(yearList, sectorNames, Emat); colormap(parula); colorbar;热力图在这里还有一个用途排错。某个行业某年数据如果异常偏大或偏小heatmap 上会出现明显色块比看数字表更容易发现问题。常见的原因包括活动数据录入单位错误、排放因子这一年被替换、以及某年缺数据被补了 0。发现异常后回头查原始表比对着公式找错快得多。3. 减排策略情景推演在MATLAB里参数化所有如果减排策略本身并不直接减它是通过改变活动水平和排放因子间接改变排放的。提高能效单位产值能耗下降活动数据变小电源结构清洁化单位电量的排放因子变小电气化率上升终端化石燃料被电力替代活动结构发生迁移。把这些策略翻译成参数模型就变成了一台可以反复按下如果键的推演器。关键在参数要少而全参数太多组合空间爆炸分析无从下手参数太少又说不清不同策略之间的差别。3.1 把减排措施翻译成可调参数的情景参数表我一般建一个 scenario 参数表而不是把参数散落在脚本里。散落变量的问题是改一处漏一处表结构让参数、取值和说明待在一起后面按行循环和输出对比都干净。基准情景取最近三年的实际均值政策情景取行业规划或可研报告里的目标值每个参数都要能说清来源。参数物理含义基准情景减排情景数据依据alpha能源强度年下降率%0.51.2行业节能目标、技改报告renew可再生能源发电占比%2035电源发展规划elecShare电力排放占总排放比例%3535由第 2 章核算结果计算beta需求侧增长不确定性%22产值预测区间elecShare 不是策略变量它是从核算模型里算出来的结构参数固定下来减小分析维度。alpha 和 renew 才是策略旋钮后面做交叉推演时只动这两个。beta 放进蒙特卡洛而不是情景表因为它描述的是不确定性而不是策略选择二者混在一起会让结果解读变得困难。3.2 用 meshgrid 展开多策略组合9 个情景一次算完两个参数各取三个水平全部组合是 9 个情景。meshgrid 是展开这种组合最直接的工具生成两个矩阵经过 I(:) 展平后拼进一张表每行一个情景。三个以上参数做交叉时组合数长得很快所以一般只对最关心的 2 到 3 个参数做交叉其余参数固定为中间值。% 两个策略参数各取 3 个水平共 9 个组合情景 alphaList [0.5, 0.8, 1.2]; % 能源强度下降率 renewList [20, 28, 35]; % 可再生能源占比 [A, R] meshgrid(alphaList, renewList); scen table(A(:), R(:), VariableNames, {alpha, renew}); % 基准年排放(万tCO2e)由第 2 章核算结果直接读入 E0 100; elecShare 0.35; Emat zeros(height(scen), 10); % 逐年推演强度下降是复利电力清洁化按占比折算 for k 1:height(scen) for t 0:9 Emat(k, t1) E0 * (1 - scen.alpha(k)/100)^t ... * (1 - scen.renew(k)/100 * elecShare); end end推演公式里(1 - alpha/100)^t是能效改进的复利效果t 从 0 到 9 表示推演 10 年。(1 - renew/100 * elecShare)表达的是可再生能源占比每提高 1 个百分点大约有 elecShare 比例的排放跟着下降。这里把电力部门和非电力部门合并成一行简化了模型代价是无法单独观察电力行业的路径需要细化时把电力排放单独拆一列而非电力部分完全不受 renew 影响。循环里 Emat 的行对应情景、列对应年份画图时直接plot(Emat)就能得到 9 条曲线组成的曲线族。3.3 蒙特卡洛抽样用拉丁超立方给减排路径画置信带点估计只能给出一条线业务方下一步一定会问如果参数没那么乐观怎么办。把 alpha、renew、beta 当成随机变量抽样跑几百次用分位数区间回答这类问题。抽样方法用 lhsdesign 而不是 rand原因是拉丁超立方把参数空间分层后再取样300 到 500 组样本就能稳定覆盖而普通随机抽样需要更多样本才能达到同样的覆盖率。n 500; % 抽样次数300~500 足够稳定 params lhsdesign(n, 3); % 3 列 [0,1] 均匀分层样本 % 把 [0,1] 映射到参数实际区间 alpha 0.5 params(:,1) * 0.9; % 能效提升率 0.5% ~ 1.4% renew 20 params(:,2) * 15; % 可再生占比 20% ~ 35% beta 0.02 params(:,3) * 0.03;% 需求侧年波动 2% ~ 5% E_sim zeros(n, 10); E0 100; elecShare 0.35; for i 1:n for t 0:9 E_sim(i, t1) E0 * (1 - alpha(i)/100)^t ... * (1 - renew(i)/100 * elecShare) * (1 beta(i))^t; end end % 逐年取 5%、50%、95% 分位数 pct prctile(E_sim, [5, 50, 95], 1);映射公式是min 样本 * (max - min)把 [0,1] 的样本平移到真实区间。prctile 的第三个参数 1 表示沿第一维样本维逐列计算得到 3×10 的结果第一行是最悲观线第二行是中位数第三行是最乐观线。画图时用 fill 填充中位数到上下界的区间带比画几十条散线清楚。另外一个细节beta 带 t 次方意味着需求侧的不确定性会随时间放大这正是蒙特卡洛要暴露出来的信息——长期预测的置信带必然比短期宽宽到什么程度是决策者该看见的。4. 组合BP神经网络与优化工具箱寻优减排策略4.1 为什么核算模型之外还要 BP 神经网络拟合曲线核算模型是白箱每一步都有物理含义但它依赖活动数据和排放因子齐全。现实中常见的情况是只有产值、人口、气温这类间接指标排放和它们的关系非单调、说不清机理这时数据驱动的黑箱拟合就有用。BP 神经网络在 MATLAB 里用 feedforwardnet 就能搭典型用途有两个一是做短期排放预测输入经济指标和气象变量输出年度排放二是做排放因子的反标定用实测序列反推等效因子。要提醒的是 BP 替代不了核算它只做拟合和预测策略比选时也有人用层次分析法AHP打分AHP 适合定权重的决策场合想讲客观最优仍然要把目标函数和约束写清楚交给优化求解更可复现。4.2 MATLAB 里训练 BP 网络的最小流程与参数设置% 输入经济产值、人口、能源消费强度输出年度排放 X [gdp; pop; intensity]; % 3×N 输入矩阵 Y emissions; % 1×N 输出向量 % 归一化到 [0,1]psX/psY 保存映射参数用于反归一化 [Xn, psX] mapminmax(X, 0, 1); [Yn, psY] mapminmax(Y, 0, 1); % 单隐藏层 10 个神经元trainlm 收敛快但吃内存 net feedforwardnet(10, trainlm); net.trainParam.epochs 1000; net.trainParam.goal 1e-5; net.trainParam.min_grad 1e-7; % 时间序列必须按顺序切分最后 20% 做测试 net.divideFcn divideblock; net.divideParam.trainRatio 0.8; net.divideParam.valRatio 0; net.divideParam.testRatio 0.2; % 训练并做全样本回代 [net, tr] train(net, Xn, Yn); Y_pred mapminmax(reverse, net(Xn), psY); test_idx tr.testInd; mape mean(abs((Y_pred(test_idx) - Y(test_idx)) ./ Y(test_idx))) * 100;输入矩阵约定每列一个样本、每行一个特征输出是行向量。mapminmax 把数据压到 [0,1] 避免神经元饱和训练完必须用 ps 结构做反归一化否则预测值和真实值根本不在一个量级。feedforwardnet(10, trainlm) 表示单隐藏层 10 个神经元trainlm 是 Levenberg-Marquardt中等样本收敛快数据噪声大时换成 trainbr 更稳。divideFcn 默认是随机切分时间序列必须改 divideblock把最后两年留作测试否则测试集混在样本中间等于提前偷看了未来信息。提示样本不到 100 个时隐藏层神经元从 5 个起步别一上来就 20。神经元越多越容易在训练集上拟合得漂亮测试集上一塌糊涂。tr 结构里存了训练和测试索引tr.best_epoch 记录最佳迭代轮次作图时可以画出训练误差和测试误差的下降曲线看两者是否在某一轮开始分叉。如果测试集 MAPE 在 5% 以下模型基本可以对外讲结论超过 10% 就要回去检查是输入变量漏了关键驱动因素还是神经元数量不对而不是急着堆数据。4.3 将减排策略寻优写成带约束的优化问题减排策略寻优的标准写法是决策变量 x_i 表示第 i 项措施的实施水平连续取值 0 到 1。目标是总成本最小约束有两个年度排放不超过上限总成本不超过预算。减排量要和核算模型的基准 E0 对应用年口径如果是跨年规划把逐年约束展开每一年的排放都写成 x 的线性函数。这类问题在 MATLAB 优化工具箱里用 optimproblem 描述最省事目标函数和约束写成数值表达式求解器会自动匹配算法。4.4 用 MATLAB 优化工具箱的 linprog 求解并检查解的合理性% 5 项减排措施实施水平 0~1 x optimvar(x, 5, LowerBound, 0, UpperBound, 1); unitCost [12; 8; 25; 40; 6]; % 单位减排成本万元/万tCO2e reduction [30; 18; 55; 90; 10]; % 各措施年减排量万tCO2e E0 100; cap 60; budget 800; % 基准排放、排放上限、预算上限 prob optimproblem(ObjectiveSense, minimize); prob.Objective unitCost * x; % 总成本最小 prob.Constraints.emissionCap E0 - reduction * x cap; prob.Constraints.budget unitCost * x budget; [sol, fval] solve(prob, Options, ... optimoptions(linprog, Display, iter));optimvar 声明决策变量并给出边界0 表示不实施1 表示满额实施。排放约束写成E0 - 减排量 cap的原因是引入措施后排放等于基准减去减排量。solve 检测到目标与约束都是线性的会自动调用 linprog如果目标或约束里出现非线性需要换成 fmincon它要求提供初值且对初值敏感建议跑多个随机起点。连续解 x 等于 0.6 表示实施到六成实际项目可以按此排期如果措施必须整体上马把 optimvar 加 Type, integer求解器会切到混合整数规划。Display, iter 打印每轮迭代信息用来判断求解过程是正常收敛还是在约束边界上反复摆动。问题特征推荐求解器需要注意的地方线性目标 线性约束linprogsolve 自动选择稳定但只能处理线性问题非线性目标或约束fmincon对初值敏感需要多起点重试混合整数/大规模组合ga 或 intlinprog规模大时耗时明显先设小种群验证求解结果先做物理性检查各项措施的水平是否都在 0 到 1 之间总成本是否恰好贴着预算或排放上限。如果某个约束明明没到边界说明这个约束不起作用可以删掉缩小问题规模如果解里某项措施是 0但它旁边有个成本更高、减排更少的选择大概率是约束写漏了回头排查约束矩阵。5. 模型验证与敏感性分析让温室气体排放模型可交付5.1 用时间序列留出法回代验证模型模型交付前必须回答凭什么信你。核算模型的验证方式是拿历史年份的活动数据回代把算出的排放总量和公开统计口径做对比误差在几个百分点内说明因子选取合理。BP 模型则靠测试集说话时间序列数据不能随机打乱切分必须按顺序留出最后两年代码里 divideFcn 设为 divideblock 就是这个目的。测试集 MAPE 低于 5% 可以对外讲结论10% 以上就是模型结构有问题继续加数据往往没用应该回去检查输入变量是否漏了关键驱动因素。5.2 局部敏感性分析找出最该优先采集的数据局部敏感性分析的做法是逐个参数上下浮动 10%观察输出变化率。归一化敏感度大于 1 说明输出被参数放大这类参数就是数据采集的优先级所在——它的误差对结果影响最大。把第 3 章的推演逻辑封装成一个函数 run_sim敏感分析、情景推演、优化寻优共用同一个模型避免优化用的是新版本、敏感性用的还是旧版本这类交付事故。params0 [0.8, 25, 0.3]; % 基准参数alpha、renew、beta delta 0.1; sens zeros(numel(params0), 1); for i 1:numel(params0) p_up params0; p_up(i) p_up(i) * (1 delta); p_dn params0; p_dn(i) p_dn(i) * (1 - delta); E_up run_sim(p_up); % 返回 10 年后排放 E_dn run_sim(p_dn); % 归一化敏感度输出相对变化 / 输入相对变化 sens(i) ((E_up - E_dn) / mean([E_up, E_dn])) / (2 * delta); end figure; barh(sens); set(gca, YTickLabel, {alpha, renew, beta}); xlabel(归一化敏感度);分母里的 2 是上下浮动造成的总相对变化这样算出的敏感度就是单位相对变化下的输出响应三个参数之间可以互相比较。alpha 和 renew 通常是最敏感的两个因为它们直接进入复利公式beta 的影响随时间放大所以推演年份越长它的敏感度排名越高。注意局部敏感性分析要求输出对参数连续。参数是离散取值的场景比如月份取整、措施只能整年上马不适合用这个方法改用全局敏感性分析逐个组合跑一遍看分布。5.3 交付前必查的四个 MATLAB 数值坑第一单位。活动数据是吨还是万吨排放因子是 tCO2/t 还是 gCO2/kWh混在一起算出来的数字能差好几个数量级。我的做法是在核算脚本开头统一一个注释块把所有单位标注清楚再写公式。第二GWP 口径。100 年和 20 年基准下 CH4 的折算值差了快三倍报告里写明了哪套口径别人才能复现你的数字。第三中文列名。readtable 不加 VariableNamingRule,preserve 时中文列名会被改成 Var1、Var2等脚本跑完才发现引用全错浪费的时间足够把代码重写一遍。第四优化初值。fmincon 换个 x0 结果可能完全不同至少跑 5 个随机起点或者先用 ga 出一组接近最优的解作为初值再交给 fmincon 精修。 交付时把敏感性分析结果和核算口径说明打包给业务方让对方标出哪些数据最没底这一步比模型精度更能建立信任。本文还有配套的精品资源点击获取