1 引言(Introduction)

非道路移动源是指空间位置随时间变化的非道路机械大气污染源(张景文等, 2017).随着我国农业机械化水平的提高, 以农用拖拉机为主的非道路移动机械保有量急剧增加(张楚莹等, 2008), 同时农用拖拉机尾气带来的大气污染问题日渐突出(Andrew et al., 2000; Zhu et al., 2011).各级政府针对该问题出台了相关法律法规并取得了一定效果, 但大气污染问题的治理仍然任重而道远.

排放清单是大气污染物排放治理的基础, 建立完善的排放清单对改善区域大气环境具有重要意义.我国移动源排放清单研究相较国外发达国家起步较晚.早在20世纪70年代, 西方发达国家就已经开始着手移动源排放模型的研发工作, 先后开发了COPERT模型、IVE模型(姚志良等, 2006)、MOBILE模型(霍红等, 2006)等主要针对道路移动源的排放模型.鉴于主流移动源排放模型均仅适用于道路移动源排放模拟, 美国环保署(United States Enviro-nmental Protection Agency, USEPA)率先开发了专门面向非道路移动源的NON-ROAD模型(林秀丽等, 2009), 并于2005年发布了功能更为全面的综合移动源排放模型MOVES(Motor Vehicle Emission Simulator)(黄冠涛等, 2010).与之相比, 国内移动源排放模型尚处在研发阶段.在排放清单的编制标准方面, 美国环保署早在20世纪90年代就提出了针对不同排放源的排放清单编制标准(U.S. Environmental Protection Agency, 2008).我国相比美国等发达国家在排放清单编制标准制定方面具有明显的后发优势, 在吸收总结国外研究经验的基础上针对移动源污染物排放特征开展了许多研究工作, 并于2014年由中国生态环境部发布了《道路机动车大气污染物排放清单编制技术指南(试行)》(http://www.mee.gov.cn/gkml/hbb/bgg/201501/W020150107594587831090.pdf)与《非道路移动源大气污染物排放清单编制技术指南(试行)》(http://www.mep.gov.cn/gkml/hbb/bgg/201501/W020150107594587960717.pdf.)(以下简称“《指南》”), 为我国移动源排放清单编制提供了参考与指导.

目前国内部分学者已经开展了基于排放清单的非道路移动源污染物排放研究, 如张礼俊等(2010)基于非道路移动源活动水平, 结合排放因子法建立了2006年珠江三角洲地区排放清单, 并分析清单数据发现非道路移动源已成为该地区第三大SO₂和NO_x排放贡献源.在长三角区域, 有学者以该地区三省一市非道路移动源实地调查数据为数据源构建排放清单, 并基于非道路移动机械保有量更新率对未来若干年的年排放总量进行预测, 为污染物排放的提前预防提供了理论依据(鲁君等, 2017; 黄成等, 2018).范武波等(2018)采用排放因子法结合非道路移动机械保有量建立了四川省排放清单, 并在此基础上对各种非道路移动机械污染物排放贡献率进行研究分析, 发现农业机械的污染物排放量仅次于工程机械, 在所有非道路移动源中排名第二.在西部地区, 姜卓等(2019)针对乌鲁木齐非道路移动源排放问题开展研究, 依据当地非道路移动机械活动水平编制排放清单, 研究发现农业机械在非道路移动机械排放量中占比最大.综上可以发现, 非道路移动源污染物排放研究对于大气污染治理具有较大意义, 但国内研究区域多集中在东部沿海, 有关中国西部地区的研究相对较少(Wang et al., 2016).

基于此, 本研究以新疆为研究区域, 采用NON-ROAD模型和《指南》提供的方法, 选取非道路移动机械中排放占比较高的农用拖拉机作为研究对象, 分别估算新疆1993—2018年PM₁₀、PM_2.5、CO、HC和NO_x等5项污染物的排放清单, 并在此基础上运用基于广义神经网络(Generalized Regression Neural Network, GRNN)与MIV(Mean Impact Value)算法构建的GRNN-MIV组合模型(吐尔逊·买买提等, 2017)定量分析排放清单时间序列的影响因素, 并结合灰色关联度分析法(Grey Relation Analysis, GRA)和斯皮尔曼等级相关系数法(Spearman Correlation Coefficient)分析污染物排放影响因素.

2 研究方法(Research method) 2.1 模型选取

目前MOVES模型是众多移动源排放模型中应用最为广泛的, 国内学者在相关领域的研究当中多次应用MOVES模型, 并取得了一定的研究成果(谢荣富等, 2017; 郝艳召等, 2017; 姚胜永等, 2019).最新发布的MOVES2014b模型针对非道路移动源增加了NON-ROAD模型(U.S. Environmental Protection Agency, 2015), 该模型内嵌大量非道路移动源排放的台架和车载实测数据, 专门应用于非道路移动源的排放模拟, 因此, 本研究选用MOVES2014b内置的NON-ROAD模型对新疆农用拖拉机的排放因子进行排放模拟与测算. NON-ROAD模型提供的非道路机械实际排放因子的具体计算公式如下:

(1)

式中, EF_adj为实际排放因子(g·hp^-1·h); EF_ss为稳定工况下排放因子(发动机0使用单位时间排放因子), 它是由实验室或生产厂家给出的新发动机在几个一般的常用的工作循环下的排放因子(g·hp^-1·h); TAF为瞬态校正系数, 表示恒定条件测定下的排放和实际运行状况下排放的比例关系; DF为劣化系数, 是指随着发动机使用时间的增加导致排放增加的系数; CF为其他的修正系数, 如温度修正、硫含量修正.具体各项修正系数计算方法参见文献(林秀丽等, 2009)

2.2 NON-ROAD参数本地化修正

MOVES模型内置的NON-ROAD模型可以根据测算的空间范围分为国家(Nation)和郡县(County), 为了使模型更加贴合研究地区的实际情况, 选择与研究区域气候条件和地理条件相似的地区作为模拟区域, 再通过Non-ROAD data importer模块实现本地相关数据的输入.本研究为减少NON-ROAD模型中默认数据对测算结果的干扰, 通过输入更符合新疆实际情况的气象参数、燃油参数、车龄分布数据与拖拉机活动水平数据, 通过NON-ROAD模型实现排放因子的本地化修正, 并在此基础上运用本地化排放因子估算污染物排放量.

时间与空间设定: 测算的起始年选择为1993年, 测算周期为1年, 根据农用拖拉机作业特点进行模拟, 依次对历年排放情况进行测算.新疆的纬度范围为34°22′~49°10′N, 地理条件较为复杂多变, 有天山山脉横亘中部, 而美国蒙大拿州(Montana)在纬度上与我国新疆较为接近, 同时也有落基山脉穿过该州, 因此, 选择蒙大拿州作为新疆的模拟地区.

拖拉机保有量和车龄分布: 拖拉机保有量数据来自1994—2019年《新疆统计年鉴》, 报废量数据来自新疆农机局并通过公式(2)核算逐年的拖拉机车龄分布.

(2)

式中, S (Y)为Y年新车数量(辆); P (Y)为Y年保有量(辆); P (Y-1)为Y-1年保有量(辆); D (Y)为Y年报废量(辆)。

拖拉机活动水平数据采用《指南》中提供参数.本研究根据新疆的实际情况, 依据《车用柴油(GB19147—2016)》中的燃油参数来进行录入.气象参数主要是温度和湿度, 气象数据来自中国气象数据网与《新疆统计年鉴》.测算污染物包括PM₁₀、PM_2.5、HC、NO_x和CO.移动源机械类型为非道路柴油农业机械.

2.3 污染物排放量计算方法

污染物排放量通过指南法与模型法提供的污染物排放因子结合新疆农用拖拉机活动水平数据计算得到.具体计算公式(http://www.mep.gov.cn/gkml/hbb/bgg/201501/W020150107594587960717.pdf; U.S. Environmental Protection Agency, 2018)如下:

(3)

式中, E为非道路移动机械各污染物排放量(t); j为非道路移动机械的类别; k为排放阶段; n为功率段; P为拖拉机保有量(辆); G为平均额定净功率(kW·台^-1); LF为负载因子; HR为年使用小时数(h); EF为污染物排放系数(g·kWh^-1).

式(3)中排放因子分别来自NON-ROAD模型模拟与《指南》推荐值, 其它参数均使用《指南》推荐参数, 农用拖拉机功率段及排放阶段分布参考文献(吐尔逊·买买提等, 2019), 拖拉机保有量数据来自《新疆统计年鉴》.

2.4 影响因子定量分析方法

GRNN-MIV组合模型是基于GRNN神经网络与MIV算法的组合模型, 本研究以该模型为主要研究方法结合灰色关联分析与斯皮尔曼等级相关系数对影响因素进行定量分析.

2.4.1 GRNN神经网络

GRNN是一种具备柔性网络结构同时兼具较强的容错性和鲁棒性的神经网络, 因此, 该神经网络很适合用于数据样本较少的非线性问题(王小川, 2013).其网络结构如图 1所示.

GRNN在结构上由4层组成, 分别为输入层(Input layer)、模式层(Pattern layer)、求和层(Summation layer)及输出层(Output layer), 在输入层输入矩阵X =[x₁, x₂, …x_n]^T, 对应输出层输出矩阵Y =[y₁, y₂, …y_k]^T.

① 本研究中影响因素样本数量较少, 因此, 选择其中16年样本作为训练集, 剩余部分作为测试集.

② GRNN神经网络中参数spread对该神经网络的性能影响较大, 选取合适的spread值对GRNN神经网络的构建意义重大.k折交叉验证法(k-fold cross-validation)是最为常用的交叉验证方法之一, 本研究利用k折交叉验证法与误差最小化方法寻求spread的最优值, 方法见文献(骆昕, 2014).

2.4.2 MIV分析法

平均影响值算法是神经网络中评价自变量对因变量影响程度最好的指标之一, MIV值的绝对值大小代表某属性影响力的大小, 正负号代表影响方向.

① 建立数据矩阵S, S是一个M行、N列的矩阵, 其中, M表示样本数量, N表示特征个数.对神经网络进行初始化操作, 将数据集S作为训练集用于训练网络.

② 网络训练完成后, 将数据集S中的每个属性分别乘以110%与90%形成新的训练集S1与S2, 并将这两个训练集输入神经网络.

③ 求两个输出结果之差, 该值即为Impact Value, 可以反映自变量对因变量的影响程度.

④ 按样本数目求得该自变量对因变量的平均MIV值.

⑤ 重复以上步骤, 分别求出每个特征的MIV值, 根据MIV值的大小对各个特征进行排序, 判断每个特征对输出结果影响的大小, 方便后续研究的开展.

2.4.3 灰色关联度分析

灰色关联分析是一种对系统发展态势进行定量描述和比较的方法, 该方法常被用于分析各个属性对于结果的相关程度(梅振国, 1992), 其值越大表示相关性越强.基本计算步骤如下:

① 确定分析数列, 影响系统行为的属性数据集为比较数列, 即子序列; 反映系统行为的数据集为参考数列, 即母序列.参考数列为: Y=Y(k), k=1, 2, 3…n; 比较数列为: X_i=X_i(k), k=1, 2, 3…n, i=1, 2, 3…m.

② 由于各个属性数据数值的数量级不同可能会使结果产生偏差, 因此, 要对数据进行无量纲化处理, 基本采用初值化或均值化的方法进行处理.

③ 计算关联系数, 计算公式如下:

(4)

(5)

式中, k为对应时间段, 即年份; i为特征; ρ为分辨系数, 通常在区间(0, 1)中取值, 本研究中ρ取0.5.

④ 分别求各属性关联系数的平均值作为子序列与母序列间关联程度的数量表示.

2.4.4 斯皮尔曼等级相关系数

斯皮尔曼等级相关可以用来估计两个变量X与Y之间的相关性, 同时, 斯皮尔曼相关性系数对数据条件要求较小, 适用性较广.其基本原理为: 对两个变量X、Y进行排序, 记录排序后X、Y中元素的位置, 记为X′与Y′, 也称之为秩次, 两者之差记为d_i, 通过d_i计算得到相关性大小.计算公式如下:

(6)

式中, N为样本数量, d_i为秩差.

3 结果与分析(Results and analysis) 3.1 新疆农田作业拖拉机污染物排放清单

利用新疆各年度拖拉机保有量结合模型法排放因子计算新疆1993—2018年农田作业拖拉机污染物排放清单, 结果见表 1(考虑到PM₁₀中含PM_2.5, 因此, 表 1总量为PM₁₀、CO、HC、NO_x之和).

3.2 排放特征分析 3.2.1 模型法排放特征及其演变趋势

表 1结果表明, 1993—2018年新疆农田作业拖拉机5种污染物年均排放量分别为0.19×10⁴、0.20×10⁴、0.19×10⁴、1.80×10⁴及1.03×10⁴ t, 26年总排放量增加了13907.85 t, 年平均增长率为2.3%.研究时域污染物排放演变曲线揭示, 1993—2018年污染物排放分为5个阶段: 第1阶段(1993—2005年)为污染物排放平稳上升阶段, 该时段新疆社会经济水平稳步提升, 农用机械普及率平稳增长, 根据新疆统计局数据, 新疆农用拖拉机数量从172322辆增加到406900辆, 增加了234578辆, 从而造成拖拉机污染物排放量的逐年提高; 第2阶段(2006—2008年)为污染物排放波动阶段, 这3年间污染物排放量出现了一定程度的震荡, 可能的原因是农机补贴、农田生产作业固有的季节性、退耕还林等可能会导致污染物排放量的波动; 第3阶段(2009—2013年)为污染物排放快速增长阶段, 该时期污染物排放量保持较高的环比增速增长, 年均增长率约为7.12%, 由于中央惠农政策的落实, 新疆第一产业结构调整加快, 农业生产由原先的劳动密集型生产模式逐步转向技术密集型的机械化生产, 使得该阶段农业机械数量增长较快, 年平均增长数量达到29748.25辆, 导致农业机械污染物排放快速增长; 第4阶段(2014—2016年)为污染物排放平稳阶段, 由于前期农业机械保有量快速增长, 农业机械数量逐渐趋于饱和, 自2014年起, 新疆农用拖拉机数量增长开始放缓, 根据新疆统计年鉴计算发现, 2014—2016年拖拉机保有量年增长率仅为1.2%, 农用拖拉机保有量增长的减缓, 导致该阶段污染物排放总量趋于平稳; 第5阶段(2017—2018年)为污染物排放下降阶段, 由于2017年是国家《大气污染防治行动计划》的收官之年, 新疆政府在《新疆维吾尔自治区2017年度大气污染防治实施计划》中明确提出2017年底前全区基本淘汰黄标车, 因此, 拖拉机数量上虽有少量上升, 但农机总体排放水平有所下降, 故而2018年农用拖拉机污染物排放量相比2017年下降了11688.79 t。但由于新疆较大的农用机械基数, 导致新疆农业机械污染物排放总体趋势依然严峻.

3.2.2 单位功率拖拉机排放趋势

单位功率排放量可以用于解释在农机总动力不断增长的背景下, 农机总动力与排放量之间的关系, 并且在排放量影响因素的研究当中也具有一定意义.其计算方法为: 单位功率排放量=当年排放总量/当年拖拉机总动力.根据表 1可知, 1993—2018年单位功率排放总体呈下降趋势, 整体下降了67%, 下降幅度较大.

单位功率排放的变化趋势说明政府层面出台的大气污染物防治措施以及基于更高排放控制技术的农业机械污染物排放标准对控制新疆农用拖拉机污染物排放具有一定的作用.

3.2.3 不同类别拖拉机排放分析

不同类型农用拖拉机在污染物排放情况与演变趋势上必然存在差异.根据图 2可知, 小型拖拉机排放量在2003年达到峰值, 随后开始呈现逐年下降的趋势, 仅在2018年出现了小幅度反弹(图 2a).与小型拖拉机排放变化趋势相比, 大中型拖拉机排放量上升幅度明显较大, 大中型拖拉机的污染物排放量逐年上升(图 2b). 2006年是两类拖拉机污染物排放量发生变化的重要时间节点, 在此之前小型拖拉机污染物排放量均高于大中型拖拉机, 两类拖拉机平均污染物排放量占比分别为59.7%与40.3%, 小型拖拉机排放量比大中型拖拉机高出近20%, 而2006年之后大中型拖拉机排放量开始反超小型拖拉机, 大中型拖拉机的年均排放量占年平均排放总量的70.9%, 而小型拖拉机年均排放量占比则为29.1%.由此可以说明新疆农业机械正在逐步实现大型化, 随着大中型拖拉机保有量不断增加, 其污染物排放占比也逐步增高.因此, 控制大中型拖拉机污染物排放将成为拖拉机污染物防控工作的重点.

此外, 分析图 2中不同污染物的分担率发现不同类型拖拉机的污染物排放分担率基本一致, 大中型拖拉机NO_x、CO、PM₁₀、HC平均分担率分别为55.1%、32.5%、6.3%、6.1%, 与之相对应的小型拖拉机分担率为5.1%、32.5%、6.2%、6.1%.从中可以看出, NO_x与CO已经成为农用拖拉机尾气排放物中最主要的污染物, 在今后大气污染防治工作中需重点关注.

3.3 指南法与模型法比较

本研究使用《指南》推荐排放因子和本地化的NON-ROAD排放因子对新疆1993—2018年的拖拉机排放进行估算.分析图 3计算结果可以发现, 基于2种方法测算的排放总量在总体变化趋势上基本相似, 但在污染物排放量上两者存在差异.1993—2003年模型法污染物排放总量计算结果与指南法相比平均高出13.7%, 但2004年起两者差异开始逐步扩大, 指南法计算结果相比模型法高出8.6%~47.9%.其中, 1993—2003年, 模型法计算所得的可吸入颗粒物(PM₁₀、PM_2.5)、HC、NO_x、CO等污染物的排放量相比指南法计算结果分别平均高出51.5%、14.6%、5.9%、14.4%.而从2004年开始指南法计算结果相比模型法则分别平均高出36%、47.7%、30.6%、35.6%.从中可以发现, 计算所得污染物排放总量时间序列的变化趋势具有一定可信度, 但两种方法的污染物排放总量计算结果存在一定差异, 其中, 2004—2018年使用《指南》推荐排放因子计算的拖拉机排放明显高于采用NON-ROAD模型本地化排放因子得到的结果.

经过分析发现造成这种差异的原因是多方面的, 而最主要的原因一方面是由于《指南》提供的排放因子是全国非道路移动源排放因子的平均水平, 并不一定完全符合新疆的实际情况, 并且《指南》编制于2014年, 将其应用于近20年前的污染物排放量计算显然缺乏时效性, 因此, 在计算新疆拖拉机排放清单时必然会出现一定的不确定性.另一方面, 在利用NON-ROAD模型本地化排放因子的过程中由于数据可获得性等因素的影响, 使得排放因子本地化不够完全, 从而导致计算结果可能与实际排放情况存在差异.

3.4 排放清单不确定性

本研究中估算的排放清单不确定性主要包括活动水平和排放因子这2个方面.①国内目前缺乏农业机械活动水平的相关文献, 因此, 本文主要参考《指南》中提供的数据, 其是基于较为宏观层面的调查统计获得, 所以并不能很好反映新疆等边缘区域的实际情况, 从而造成了一定的不确定性.②在利用NON-ROAD模型对排放因子本地化的过程中本地化处理仍不够全面, 由于数据可获得性等因素的限制, 排放因子计算过程中较为关键的稳态工况排放因子仍采用了NON-ROAD模型提供的几种机型拖拉机的排放因子, 本研究对这种由国内外机型差别对排放因子带来的影响未做深入考虑。

3.5 拖拉机排放影响因素分析 3.5.1 影响因素选取

非道路移动源污染物排放数据是具有季节性、混沌性和随机性特征的时间序列, 因此分析其影响因素时, 本文从农业经济发展水平、农业生产作业水平2个维度出发, 结合文献(张乐勤等, 2012; 陈亮等, 2018)和专家咨询法建立由农机总动力(T1)、大中型拖拉机保有量(T2)、小型拖拉机保有量(T3)、拖拉机年增长量(T4)、农作物机耕面积(T5)、农业生产总值(T6)、农村个人固定资产投资(T7)及污染物排放总量等8个指标组成的污染物排放影响因素关系数据模型.T1~T7均数据来自《新疆统计年鉴》.

3.5.2 影响因素量化

本研究采用灰色关联度分析法(GRA)与斯皮尔曼等级法(Spearman)(具体方法与步骤见2.4.3节与2.4.4节)分别将各影响因素与污染物排放总量的相关程度进行量化(表 2), 并与GRNN-MIV模型结果相结合进行影响因素分析.

为进一步探究各影响因素对污染物排放总量的影响程度, 本研究将T1~T7作为输入, 污染物排放总量作为输出, 将2000—2016年的数据集作为训练集, 剩余作为测试集, 基于Matlab平台编程实现GRNN-MIV模型.神经网络实验结果表明, 当spread=0.56时, 验证集预测误差最小, 因此, 以spread值为0.56时的网络输入和输出为输入输出神经元构建最终GRNN网络, 并选用MSE、MAPE和R-Square 3项指标作为神经网络的检验指标, 对2000—2016年污染物排放总量进行拟合.验证结果表明, 拟合的MSE、MAPE和R-Square分别为164.73、1.39%和0.99.从3项检验参数的大小可以发现, 模型拟合准度较高, 可以用于污染物排放总量影响因素分析.

基于上述神经网络, 以2000—2018年的影响因素数据与排放总量数据作为数据源, 计算各影响因素的MIV值, 结果如图 4所示.采用标准化方法对MIV值进行处理, 标准化后的值分别为0.41、0.15、0.13、0、1、0.25、0.09.

结合表 2与图 4分析发现: ①农作物机耕面积对农用拖拉机污染物排放总量的影响最为显著, 3种方法的计算结果均显示农作物机耕面积与污染物排放总量具有强相关性. ②农机总动力、农业生产总值、农村个人固定资产投资作为农业经济发展水平的重要指标MIV值分别为0.41、0.25、0.09, 相关性分析结果则均高于0.8, 这表明上述3项指标对污染物排放量具有一定影响.③相比小型拖拉机保有量, 大中型拖拉机保有量对污染物排放总量的影响程度更高, 这也进一步说明控制大中型拖拉机污染物排放是解决农用拖拉机污染物排放问题的关键.

4 结论(Conclusions)

1) 本文采用指南法与模型法分别计算新疆农业生产用拖拉机排放PM₁₀、PM_2.5、HC、NO_x和CO等5种污染物的清单, 结合新疆社会经济现状对比2种方法适用性发现, 基于本地化修正的NON-ROAD模型在估算研究区域污染物排放方面具有一定参考意义.

2) 1993—2018年农用拖拉机排放量总体上呈现出波动中增长趋势, 其中, 1993—2005年是污染物排放量的平稳增长期, 污染物排放量逐年缓慢增长; 2006—2008年3年间可能受到多种外界因素影响, 导致该时间段污染物排放总量出现了一定程度的震荡; 2009—2013年是污染物排放量的快速增长期, 年均增长率达到7.12%;2014—2016年污染物排放量仅有小幅波动, 总体态势平稳; 2017—2018年, 由于政府层面相关措施的落实, 污染物排放量出现了一定幅度的下降.但由于新疆农用机械保有量较大, 农用机械污染物排放总量还处在一个较高水平, 因此, 新疆的大气污染治理工作仍然任重而道远.

3) 对比26年大中型拖拉机与小型拖拉机的污染物排放分担率, 结果表明, 2006年之后大中型拖拉机年均排放量超过小型拖拉机年均排放量, 大型拖拉机年均排放量占年平均排放总量的70.9%, 小型拖拉机年均排放量占29.1%.农用大中型拖拉机是农用拖拉机污染物排放的主要贡献者, 因此, 减少大中型拖拉机污染物排放需要结合农业生产需求, 在技术层面和排放治理措施方面采取更有力举措.

4) 基于MIV算法、灰色关联度分析法和斯皮尔曼等级相关系数对污染物排放影响因素进行量化, 结果表明, 农业经济发展水平对区域污染物排放总量及变化起决定性作用, 提出的3种影响因素量化方法的组合在影响因素量化方面具有一定参考价值.

本次研究虽通过排放因子本地化降低了排放清单的不确定性, 但受数据可获得性限制, 所建立的排放清单中仍存在一些目前尚难以消除的不确定性, 在后续工作中, 将在实际工况下的排放因子测度和特性分析方面做进一步研究.