基于MATLAB的多区域投入产出模型碳排放核算实战指南

简介基于多区域投入产出模型MRIO的二氧化碳排放量计算是环境经济学与碳足迹研究的重要方法。这份MATLAB代码资源面向需要核算全球或区域间碳排放的研究者、学生及政策分析人员围绕WIODr16数据库实现了从数据预处理、Leontief逆矩阵运算到生产者责任PBA与消费者责任CBA核算的完整流程。压缩包共2个文件以Matlab的.m脚本为核心代码附带一个txt说明文件整体仅3KB轻量易用。通过这份代码使用者可以快速掌握MRIO框架下CO2排放的计算逻辑复现并扩展具体的碳足迹核算过程为减排政策或供应链绿色转型提供量化支持。目前该资源已有1194人学习适合具备一定投入产出分析基础、希望借助MATLAB开展环境核算实践的相关人群。 做环境经济核算的同行提到多区域投入产出模型MRIO应该不陌生。我前阵子刚完成一个用MATLAB实现的多区域碳排放核算任务把Eora和WIOD两套数据都跑了一遍踩了不少坑。这个项目的核心痛点很典型单独看某个区域的直接排放很简单但现代经济体系里产品和服务在区域间频繁流动你用的每一度电、买的每一件商品背后牵动的排放都在多个区域之间转移。MRIO模型正好能把这笔看不见的碳账算清楚而MATLAB凭借矩阵运算和数据处理优势是跑这套模型最顺手的工具之一。这套内容适合三类人第一类是写论文的硕士博士尤其是环境经济学、区域经济、产业生态学方向第二类是做碳核算、碳足迹、双碳咨询的从业者需要算清楚区域间贸易隐含碳第三类是单纯想学MATLAB处理多区域矩阵数据、熟悉csv导入和table操作的同学。下面我把完整思路、数据预处理方法、核心计算流程和排错经验一次讲透。1. 核心思路为什么用多区域投入产出模型算碳排放1.1 投入产出模型的基本原理投入产出模型最早由列昂惕夫提出核心思想是记录经济系统中各产业部门之间的投入与产出关系。一张投入产出表里行方向表示某个部门的产品分配给哪些部门做中间投入、哪些做最终使用列方向表示某个部门在生产过程中消耗了哪些部门的中间产品以及增加了多少劳动、资本等初始投入。以单区域模型为例总产出向量x、中间需求矩阵Z、最终需求向量y之间满足基本的行平衡关系x Z × 1 y这里的1是单位列向量Z × 1表示各行业中间使用的合计。接下来定义直接消耗系数矩阵A元素a_ij表示第j部门生产单位产品需要直接消耗第i部门的产品数量A Z × diag(1./x)然后利用列昂惕夫逆矩阵B (I - A)^(-1)可以得到完全需求系数。这也是整个投入产出分析最核心的公式它把直接消耗和间接消耗全部纳入考量。比如汽车制造要直接消耗钢铁钢铁又要消耗铁矿石铁矿石开采还要用电这一串连锁反应都能通过逆矩阵完整刻画出来。正是靠这个传导机制投入产出模型才能算清一个产品在全生命链条上的资源消耗和排放。1.2 单区域转多区域MRIO怎么处理区域间贸易单区域模型只描述一个经济体内部的部门关系但实际生产中产品中间品经常跨区域流动。比如A区域的电子元件运到B区域组装成整机整机又出口到C区域消费。这种隐含在贸易中的碳排放转移单区域模型完全无法刻画。多区域投入产出模型MRIO的做法是把各个区域的投入产出表通过区域间贸易流链接起来。假设有m个区域、n个部门那么MRIO的Z矩阵维度就是(mn) × (mn)每个块Z_sr表示s区域产品被r区域生产部门消耗的数量。同样最终需求矩阵Y的维度是(m*n) × m每一列对应一个区域的最终需求。MRIO把单区域的I-A结构扩展为多区域块矩阵对角线上的块A_ss是s区域自身的直接消耗系数非对角线块A_sr则刻画s区域产品对r区域生产部门的中间投入。这样做的好处是能区分生产中体现的排放到底是由哪个区域的最终需求拉动的也就是经常说的消费侧核算和生产侧核算。用公式表达就是完全碳排放强度矩阵F f × (I - A)^(-1)其中f是各区域分部门的直接碳排放强度向量F中的元素表示某区域某部门每生产单位最终产品在整个多区域生产网络中引起的直接和间接碳排放总量。这个思路虽然计算量比单区域大不少但在MATLAB里也就是几个矩阵运算指令的事。2. 数据准备与预处理最费时间的环节2.1 主流MRIO数据库的选择与对比MRIO计算依赖高质量的多区域投入产出表和配套的环境账户数据。目前常用数据库有Eora、WIOD、GTAP和区域级的CEADs等核心差异如下表所示。数据库覆盖范围部门个数时间跨度环境账户适用场景Eora全球约190个国家/地区26个部门含服务细分1990-2022完整含CO2等多种排放全球尺度、国家间贸易隐含碳WIOD全球43个国家/地区其他区域56个部门ISIC Rev.42000-2014含CO2、CH4等不同版本口径有差异欧盟及主要经济体、制造业分析GTAP全球140国家/地区57个部门各基准年需额外配套排放数据CGE模型衔接、贸易政策模拟CEADs中国省级尺度42部门左右2000-2019省级排放清单省际贸易隐含碳、区域减排责任如果你做全球尺度的分析优先选Eora或WIOD。Eora优势是时间序列长、覆盖全缺点是数据更新频繁且格式复杂需要仔细读取官方说明。WIOD数据口径比较规范但时间更新相对滞后。如果做中国省际研究CEADs是主流选择它的MRIO表和中国投入产出表部门口径高度一致。2.2 环境账户数据的单位与口径对齐数据下载回来后最头疼的是单位不一致。Eora的环境账户中CO2排放一般以吨为单位而经济量的单位可能因国家不同而不同有的按本币有的统一按美元。WIOD的环境账户中CO2单位是千吨经济量按百万美元计。GTAP数据库本身不含可直接使用的排放数据需要自己从EDGAR、CDIAC等排放清单中整理补充工作量更大。我自己的处理习惯是先在Excel里做一遍初步检查确认以下几点再进MATLAB所有区域的经济单位统一最好都折算成同一货币和同一价格基准。排放数据单位统一通常换算成万吨CO2或百万吨CO2避免数字过大影响后续绘图。部门口径对齐。这是最麻烦的一步比如Eora的26个部门分类和GTAP的57部门分类完全不同若做跨库对比需要先按双方部门说明做映射合并。确认MRIO表的行总计等于列总计关系是否闭合。如果数据源自带检验脚本务必先跑一遍原始校验否则后面算出来的完全消耗系数可能是错的。2.3 MATLAB导入不同格式数据多数MRIO数据库提供CSV或Excel格式MATLAB读取建议优先用readmatrix或readtable两者对混合数据类型兼容性较好。示例代码如下% 读取区域间交易矩阵Z矩阵 Z readmatrix(Eora_MRIO_Z.csv); % 读取最终需求矩阵Y Y readmatrix(Eora_MRIO_Y.csv); % 读取各区域各行业总产出向量x x readmatrix(Eora_MRIO_x.csv); x x(:); % 读取排放向量记为e e readmatrix(Eora_MRIO_CO2.csv); e e(:); % 检查维度是否匹配 disp(size(Z)); disp(size(Y)); disp(length(x));这里有个细节部分CSV文件第一列是区域或部门名称readmatrix会自动忽略非数值内容但如果有空单元格或字符串混杂最好先用readtable观察一下再用table2array提取数值部分。遇到中文表头导致乱码时使用readtable时指定编码格式能解决大部分问题。3. MATLAB核心计算流程与关键代码3.1 从Z矩阵到直接消耗系数和列昂惕夫逆矩阵数据导入后第一步是清洗维度。假设总区域数m190Eora部门数n26那么Z矩阵应为(19026) × (19026)也就是4940 × 4940。这个规模在MATLAB里直接做矩阵运算没有任何压力但如果用GTAP这类300多个区域、57个部门的超大规模数据就要考虑稀疏矩阵和分块处理。直接消耗系数A的计算要特别注意分母是各部门总产出不是总产值。MATLAB代码% 每列对应的总产出 x_col x(:); % 计算直接消耗系数矩阵A用列除 A Z ./ x_col; % 自动广播MATLAB R2016b以后支持 % 检查A中是否出现NaN或Inf fprintf(NaN数量: %d\n, sum(isnan(A(:)))); fprintf(Inf数量: %d\n, sum(isinf(A(:))));这里要强调Z中某些格子可能是0导致A中为0这是正常的。但如果x_col中出现0那对应列会出现NaN或Inf说明原始数据中该部门的产出缺失需要回到数据源核对绝对不能直接跳过。接下来求列昂惕夫逆矩阵n_total size(Z, 1); I eye(n_total); B inv(I - A); % 更推荐的写法是用左除数值稳定性更好 % B (I - A) \ I;用inv还是用左除是个细节。对中型矩阵两者差别不大但数据规模大时左除数值更稳定、速度更快。我测试过4940 × 4940的稠密矩阵左除法求逆耗时约十几秒完全可接受。如果你处理的是区域间高维MRIO可以考虑用稀疏矩阵存储I-A以降低内存占用。3.2 直接排放强度、完全排放强度和三类核算口径有了逆矩阵B就可以算三条主线。第一直接排放强度f。定义是各区域各行业每单位产出直接排放的CO2量f e ./ x;注意e和x的单位要保持一致。如果e是吨CO2x是万美元那f的单位就是吨CO2/万美元。第二完全碳排放强度F。公式为F f × B即F f * B; % 维度 1 × n_totalF中每个元素代表该部门提供一单位最终产品时整个多区域生产网络中的直接和间接碳排放总和。这个指标在做碳足迹分析时非常常见。第三把F和最终需求矩阵Y相乘能得到各区域最终需求拉动的碳排放% Yn_total × m 的最终需求矩阵 E_demand F * Y; % 维度 1 × m即各区域最终需求引起的碳排放这是消费侧碳排放核算的基础。如果想做生产侧碳排放直接对f和各部门产出求和即可E_production sum(e); % 所有区域所有部门的直接排放总和更进一步的区域间贸易隐含碳分解则要拆开块矩阵。比如从区域s流入区域r的产品所隐含的碳排放可以用以下思路计算% 拆块假设每个区域部门数相同n区域总数m for s 1:m for r 1:m rows_s (s-1)*n 1 : s*n; cols_r (r-1)*n 1 : r*n; Z_sr Z(rows_s, cols_r); % 结合A和B的块结构计算s到r的隐含碳 % 具体公式依研究口径而定 end end这个双循环在区域数多时效率较低如果只是做整体核算不推荐展开到区域对级别但若是研究贸易隐含碳流向这一步不可省。MATLAB的矩阵分块操作和parfor并行循环可以显著加速值得尝试。3.3 完整核算脚本一个可直接修改的最小示例以下是一个完整的最小化核算流程我以三区域两部门简表为例展示实际数据只需按对应变量替换即可。% 参数设置 m 3; % 区域数 n 2; % 部门数 % 模拟数据实际使用替换为读入数据 Z [10 20 30 10 5 15; ... 8 12 20 15 10 10; ... 5 10 15 10 8 12; ... 6 8 10 12 9 14; ... 4 6 8 10 7 9; ... 2 4 6 8 5 7]; Y [30 10 5; ... 20 15 10; ... 10 20 8; ... 12 8 15; ... 6 10 12; ... 4 6 10]; x sum(Z, 2) sum(Y, 2); % 总产出 中间需求 最终需求 e [50; 40; 30; 20; 15; 10]; % 各部门CO2排放单位万吨 % 计算直接消耗系数 A Z ./ x; % 计算逆矩阵 B inv(eye(m*n) - A); % 直接排放强度 f e ./ x; % 完全排放强度 F f * B; % 各区域最终需求拉动的碳排放 E_demand F * Y; % 显示结果 disp(各区域最终需求拉动的碳排放万吨CO2:); disp(E_demand);这套代码逻辑很简单但足够跑通整条核算链路。在实际项目中你可以把读数据部分替换为真实CSV路径再在末尾加上导出结果和绘图代码即可。4. 结果可视化与结果校验4.1 用MATLAB画区域碳排放结构图核算完成后至少要输出两类图一是各区域消费侧碳排放对比柱状图二是各区域各行业的完全碳排放强度热力图。柱状图用bar即可热力图用heatmap可以快速生成区域×行业的碳排放矩阵。% 柱状图各区域消费侧碳排放 figure; bar(E_demand, FaceColor, [0.2 0.6 0.8]); xlabel(区域编号); ylabel(碳排放量万吨CO2); title(各区域最终需求拉动的碳排放(消费侧)); grid on; % 热力图完全碳排放强度矩阵m*n 重排为 m×n F_mat reshape(F, n, m); figure; heatmap(F_mat); xlabel(部门); ylabel(区域); title(完全碳排放强度吨CO2/万美元);绘图完成后建议用exportgraphics直接输出高清PNG或PDFMATLAB R2020a之后的版本对矢量图导出支持已经比较好论文投稿时能用上。4.2 结果合理性校验的方法计算过程跑通了不代表结果正确我每次都会做几项核验。第一是全局平衡测试。把各部门总需求向量代入模型看计算出的总产出是否等于原始总产出x_calc B * sum(Y, 2); % 按最终需求推算总产出 diff_x max(abs(x_calc - x)); fprintf(总产出最大偏差: %.6e\n, diff_x);如果这个偏差大于1e-6说明Z矩阵或Y矩阵的数据关系有问题常见原因是原始数据未做平衡调整。第二是估算合理区间。以中国省级研究为例各省消费侧碳排放不应为负数数值量级也应落在可解释范围内。如果出现某区域消费侧碳排放远大于全球总排放就要检查f和B的维度是否搞反了。第三是跟公开结果交叉验证。比如Eora数据库官网自带一些汇总结果可以把自己计算的国家总排放跟官方发布数据做对比差距在5%-10%以内基本可接受。若偏差很大优先怀疑单位换算出错或年份选择不一致。5. 常见问题与排查技巧实录5.1 维度不匹配与CSV读取乱码MRIO数据处理中最常见的就是维度不对。下载的Z矩阵可能是m行n列也可能第一行第一列是区域部门标签。我的经验是先用size查维度再抽样检查几个已知数值是否和官方说明一致确认无误后再进行计算。避免等算完逆矩阵才发现方向反了那会浪费很多时间。CSV乱码问题在读取中文表头时经常出现MATLAB中可用opts detectImportOptions(file.csv); opts.DataLines 2; % 跳过表头 T readtable(file.csv, opts);或者直接对文件名编码处理。用readmatrix时如果遇到字符型缺失值它会自动转成NaN这一点要注意因为后续运算中NaN会污染结果。5.2 单位换算和负值处理单位换算属于低级错误但影响巨大。Eora中一个区域比如美国的CO2排放可能显示为5410百万吨而另一个数据库中是5.41×10^9吨如果没统一单位就代入计算结果会差三个数量级。另外MRIO表中可能因为统计口径不同出现负值尤其是库存变化和进口列。处理原则是负库存需要保留因为它是经济系统的一部分但对某些负值极大的格子要做逻辑判断防止A矩阵出现异常负系数。简单处理方式是把异常值设为0或按比例分摊但要在论文或报告里注明处理方式。5.3 矩阵求逆慢、内存溢出和license相关环境问题中等规模MRIO直接用inv没问题但遇到GTAP这种大规模数据内存占用可能达到几十GB。一个折中方案是使用迭代法或分块求逆% 使用稀疏矩阵优化 I_A_sparse sparse(I - A); B_sparse I_A_sparse \ speye(size(I_A_sparse, 1));此外我发现部分MATLAB版本在远程桌面环境或未激活工具箱的情况下会出现矩阵运算报错遇到这种情况优先检查系统环境和工具箱license是否正常。MATLAB下载安装后确认好license和路径设置再开始跑项目不然后续调试的时候容易分心。6. 最后分享几点实操心得跑完这套流程之后我有几个比较深的体会。第一数据预处理的时间至少要留出整个项目的60%真正的矩阵计算在MATLAB里反而不是瓶颈。Eora数据从下载到清洗成可直接计算的Excel表我花了差不多两个整天其中大部分时间花在区域代码映射和部门对齐上。建议提前画好数据流图把每一步的输入输出维度记清楚会大大减少返工。第二反向验证非常值得做。我第一次算完全排放强度后总觉得个别国家的数字偏大后来发现是环境账户的排放量单位少乘了1000。只要把计算结果跟官方发布数据对上后续所有分析才站得住脚。如果发现某区域排放结果跟常识差距太大优先检查数据而不是急着调整模型。第三如果后续还想做分解分析如SDA结构分解或LMDI指数分解建议从一开始就把区域代码和部门代码作为单独变量保存下来不要只存矩阵。因为分解分析需要频繁按区域和部门筛选数据没有明确标签会非常痛苦。这套基于多区域投入产出模型和MATLAB的碳排放核算流程整体就是数据准备、矩阵计算、结果校验三板斧。把基础链路跑通之后再往里面加行业拆分、政策情景模拟或者更细的贸易隐含碳流向分析都会顺手很多。本文还有配套的精品资源点击获取