中国田块尺度露天生物质燃烧NOx遥感排放清单

沈永林,江昌民,肖泽敏,姚凌,秦凯

(1中国地质大学(武汉)国家地理信息系统工程技术研究中心,湖北 武汉 430074;2中国地质大学(武汉)地理与信息工程学院,湖北 武汉 430074;3中国科学院地理科学与资源研究所资源与环境信息系统国家重点实验室,北京 100101;4中国矿业大学环境与测绘学院,江苏 徐州 221116)

0 引言

氮氧化物(NOx≡NO+NO2)作为对流层臭氧、羟基自由基等主要大气氧化剂以及硝酸盐气溶胶的前体物,以多种物理化学形式存在于大气中,是衡量对流层氧化能力的重要指标,也是硝酸型酸雨、光化学烟雾的主要诱因[1]。生物质燃烧是NOx排放的重要来源,研制精准的排放清单,对理解对流层臭氧、二次气溶胶的形成、氮的生物地球化学循环、农业生产合理规划布局、生态系统恢复等都具有重要意义[2-6]。

目前有多种计算区域尺度秸秆焚烧的大气污染物排放量的方法,传统生物质燃烧NOx排放估算主要采用自下而上的方法,通过地面统计资料估算生物质消耗,从而计算得到污染物排放量[7-10]。此类方法计算过程简单,能直观地反映污染物的排放水平。但排放量多根据粗粒度土地覆盖类型产品,空间分辨率低且分类尺度不统一,无法准确地反映排放源的空间分布。遥感技术具有实时、覆盖范围广的特点,时间和空间尺度上的优势使其在生物质燃烧排放估算中被广泛运用[11]。例如,Selier和Crutzen[12]提出了基于火烧迹地面积(BA)的生物质燃烧排放估算方法,并且目前已存在诸多BA的公共产品,如全球烧毁区卫星产品(GBA2000)[13]、中分辨率成像光谱仪(MODIS)BA产品[14]、全球火情数据库(GFED)[15]等。但在高纬度地区,可见光遥感对于火烧迹地的判别常受到诸如植被冠层高度、太阳高度角等因素的影响,估算误差明显,并且小型火灾极易被忽略。Kufman等[16]认为在热红外遥感中生物质燃烧产生的热辐射功率(FRP)与生物量消耗速率成比例。在此基础上,Vermote等[17]提出了基于FRP的生物质燃烧排放估算方法,避免了由小型火灾燃烧面积难以监测引起的不确定性。此后,Ellicott等[18]基于MODIS FRP观测数据估算了全球2001–2007年的火点辐射能量(FRE),并利用高斯函数模拟了全球不同区域生物质燃烧时的FRE变化。Wu等[19]利用该方法建立了我国中东部地区2003–2015年的生物质燃烧NOx排放清单;Rafee等[20]估算了北美地区的NOx排放,并模拟了其在大气传输作用下的空间分布和对环境的影响。而FRP方法的不确定性主要来源于传感器对火点光谱特征的敏感度、燃烧效率参数、排放因子参数等,故Freeborn等[21]通过连续观测位于星下点的火点,分别测算了孤立火点像元与火点像元簇FRP的标准差和均值,结果证明了该方法的可靠性。然而,当前排放清单的空间特征尺度大多停留在粗粒度土地覆盖类型燃烧引起的差异,少有考虑到精细的作物类别及其空间分布特征的重要性。

考虑到我国农业种植结构的复杂分布,本文提出了一种顾及火点像元地物类别先验概率的排放因子加权策略,并结合FRP昼夜周期分布的假设,建立中国逐日1 km的田块尺度露天生物质燃烧NOx遥感排放清单(FOBBEINOx)。即在FRE估算生物质燃烧NOx排放的框架下,利用MODIS温度异常产品(MOD/MYD14A1),结合地表覆盖类型产品(MCD12Q1)和主要农作物空间分布产品(GFSAD1000)[20]及不同作物对应的物候期,实现对我国2003–2020年生物质燃烧NOx排放估算,并将结果分别与Himawari-8数据、站点观测资料、GFED数据库以及臭氧监测仪(OMI)NO2垂直柱浓度(VCDs)进行对比验证,以评估本文方法的有效性。

1 实验方案与装置

本研究以MOD/MYD14A1温度异常产品为数据源,通过改进的高斯函数模拟地表生物质燃烧FRP的昼夜周期变化,估算生物质燃烧过程释放的能量;然后,根据我国主要作物物候期的特点,结合不同土地覆盖产品制定排放因子加权策略;最终,建立我国田块尺度露天生物质燃烧NOx排放清单(见图1)。

图1 流程图Fig.1 Flow chart

1.1 露天生物质燃烧辐射能量遥感估算方法

露天生物质燃烧FRE是焚烧时段内FRP随时间的积分。由于MODIS每日仅提供4次FRP观测,而火点燃烧最为剧烈时段的FRP常被MODIS传感器忽略,因此难以估算FRE。为解决该问题,Vermote等[17]假设生物质燃烧的FRP昼夜周期遵循高斯分布。在昼夜周期内及时整合FRP能在一定程度上弥补火灾强度被低估的情况。因此本研究根据我国农事实际情况对高斯参数进行了调整,即

式中PFR,p代表FRP昼夜周期变化的峰值,b表示背景水平,t表示时间,h代表峰值的当地时间(单位:h),σ为昼夜高斯分布曲线的宽度。Ellicott等[18]和Vermote等[17]均认为参数b、σ、h与MODIS上下午星FRP的比率有关,即Terra FRP与Aqua FRP的比(简写T/A)。Vadrevu等[22]认为不同的T/A比值能代表不同火点类型和相应FRP昼夜周期变化,并根据T/A直接计算出参数b、σ和h,具体计算公式分别为

式中x表示T/A均值。通过在公式中增加参数λ(单位:h),即可根据我国各地农业秸秆焚烧的实际情况修正昼夜曲线FRP高峰时的时间。Liu等[9]认为利用白天Aqua MODIS观测的平均PFR,Ad与b、σ、h参数,可估算出FRP峰值,其计算公式为

1.2 顾及火点像元地物类别先验概率的排放因子加权策略

排放因子是表征NOx单位排放量的关键参数。由于我国作物种植模式为庞大的、分散的小农生产,农田结构与位置复杂,田块呈现小、散、多的特点,遥感农田像元普遍存在类别混合的现象。因此,需要根据作物的空间分布、物候特性以及农业统计年鉴,设计适合中国的像元混合排放因子。本研究中各像元的排放因子由MCD12Q1和GFSAD1000共同确定。其中MCD12Q1是由国际陆界生物圈方案提供的由17种地物类型组成的500 m空间分辨率的全球土地覆盖类型方案。GFSAD1000由美国国家宇航局全球粮食供应计划提供,空间分辨率为1 km,其将超过全球耕地面积60%的作物分为8类混合作物类型,即水稻、小麦、大麦、玉米、大豆、棉花、甘蔗、木薯、果园。本研究中各地物类别的排放因子主要参考GFED v4数据库和已有研究成果[7,23-31],见表1。

表1 主要地表覆盖类型的排放因子Table 1 Emission factors of major land cover types

为有效解决像元排放因子计算过程中地物类型混合的问题,本研究设计了一种顾及火点像元地物类别先验概率的排放因子加权策略。每个像元排放因子的确定分为两步:

1)计算像元中农田的排放因子EMC。首先根据我国统计年鉴中作物类型对应省份播种面积Acrop,i(i表示作物类别)和GFSAD1000中相应作物的像元数量Pcrop,i的比值作为权重,以衡量像元中每一类作物对该混合像元亚排放因子的贡献;此外,根据作物在该地区的物候期Pcrop判断是否存在该类作物的秸秆焚烧。最后,根据GFSAD1000混合像元中所有作物类型排放因子Ecrop,i进行加权,确定作物排放因子,计算公式为

2)计算像元综合排放因子E p。当MCD12Q1将某像元区分为农田亚类,且GFSAD1000也判定该像元为农作物时,则可确定该像元为作物,即E p=EMC;当GFSAD1000判定像元为非农作物,同时MCD12Q1也为非农田,则根据MCD12Q1中该像元地物类型计算像元最终的排放因子,即E p=Elc;若GFSAD1000判定像元为作物,MCD12Q1判定为非农田,则根据MCD12Q1提供的非作物地物类型覆盖率作为权重R,而1-R作为GFSAD1000的权重,综合反映两种产品的总体分布特征。[例如:MCD12Q1产品中的某像元为阔叶林,且覆盖度大于60%,而在GFSAD1000中被认为是农作物,则MCD12Q1的权重为0.6,相应GFSAD1000的权重为0.4(见表2)]。具体的综合排放因子计算公式为

表2 排放因子权重设置规则Table 2 Emission factor weight setting rules

1.3 田块尺度生物质燃烧NOx排放估算模型

本研究根据估算的FRE和混合排放因子进行田块尺度NOx排放估算。特定地表覆盖类型生物质燃烧产生的NOx排放量Mc可表示为

式中c表示地表覆盖类型,CR为能量到质量的燃烧转化率,ENOX,c为特定地表覆盖类型的排放因子。Wooster等[32]和Freeborn等[21]分别测得CR为(0.368±0.015)kg·MJ-1和(0.453±0.068)kg·MJ-1。本研究取两者均值(0.41±0.04)kg·MJ-1,通过与FRE的乘积转换为露天生物质燃烧总质量。

2 2003–2020年中国田块尺度露天生物质燃烧NOx排放清单

根据本研究的方法可估算逐日1 km的露天生物质燃烧NOx排放量,亦可仅估算农业源。本研究统计了2003–2020年各省农业源露天生物质燃烧NOx年排放量(见表3,其中1～18分别代表2003–2020年)。广东、河北、黑龙江、吉林、内蒙古、新疆的农业源露天生物质燃烧NOx年排放量与中国的其他省份相比居高,且年际变化存在差异。2015年中共中央办公厅、国务院办公厅印发了《党政领导干部生态环境损害责任追究办法(试行)》,翌年全国农业源NOx排放量达到历年最低值。本研究排放清单结果显示这些省份政策实施效果显著,且这些地区为我国的粮食主产区[33],表明本文研究成果不仅能对秸秆焚烧政策的实施效果提供参考和依据,且能从侧面反映我国粮食主产区农业种植模式的时空变化。而河南、江苏、四川、山东等粮食主产区的农业源露天生物质燃烧NOx年排放量在2.15 Gg以下,这与其较低的排放水平与农业种植结构及对作物秸秆的处理方式有关。例如,山东处理小麦秸秆的方式除直接焚烧外,大多都将秸秆制作扫把等编制品;此外,也会将秸秆用作牛羊饲料等,且近些年来加大了秸秆禁烧工作的力度。农业源露天生物质燃烧NOx年排放量较低的还有北京、天津、上海等直辖市,主要因其城市功能定位,粮食生产功能相对较弱;同理,浙江和福建为多丘陵山区,海南、西藏、青海和宁夏粮食产量受粮食播种面积少和单产较低的双重影响。

表3 2003–2020年各省年度田块尺度农业源露天生物质燃烧NOx排放量Table 3 Annual field-level open biomass burning NOx emissions from agricultural source of each province in China during 2003–2020

Continued

图2为2003–2020年不同季节的农业源露天生物质燃烧NOx排放量。农业源NOx排放年内变化总体上呈现“U”型分布,历年在冬春两季排放量明显高于夏秋。由于我国大部分春播作物的收获期集中于9月底至11月初,当田块有播种越冬作物的需求时,如种植两至三季的麦田,秸秆会在较短的间耕期内被焚烧,以达到迅速清理耕地的目的。此外没有播种越冬作物的农田,秋收后的秸秆则可能会在田间堆积晾晒,直至第二年春播前焚烧。因此农田中秸秆焚烧季主要分布在冬春两季,夏秋季焚烧的现象则相对稀疏。

图2 2003–2020年NOx排放。(a)年排放;(b)季排放Fig.2 NOx emissions during 2003–2020.(a)Yearly;(b)quarterly

3 对比验证

3.1 与Himawari-8卫星观测的对比结果

为评估本研究中基于MODIS观测数据建立的高斯模型与实际地表生物质燃烧FRP昼夜变化的拟合优度,首先利用高时间分辨率的Himawari-8亮度温度产品数据反演得到观测频率为15 min的FRP结果[34,35],然后将反演结果与本研究模拟结果进行比较。根据MOD/MYD14A1数据中火点的地理坐标,选取湖北省境内的3例生物质燃烧实例:2016年2月7日(29.55 N,115.05 E)、2月18日(29.75 N,114.85 E)以及2月27日(29.15 N,112.95 E)。如图3所示,根据农业源生物质燃烧FRP实际分布特征,FRP峰值处于16:00–20:00,最大FRP约为1500 Mw·km2。

图3 基于MODIS观测的2016年FRP高斯分布与Himawari-8反演结果的比较。(a)2月7日;(b)2月18日;(c)2月27日Fig.3 Evaluation of 2016 FRP between Himawari-8 and MODIS observations with assumed diurnal cycle.(a)February 7;(b)February 18;(c)February 27

此外,本研究选取了2016年2月Himawari-8在湖北省及其周边的274次火点观测为样本(包括同一火点的连续观测),分别计算了对应基于MODIS模拟结果的FRP,并建立了两组数据的散点图,其皮尔逊相关系数为0.65(图4)。结果显示本研究依据MODIS昼夜FRP观测数据建立的高斯分布一致性较好。

图4 2016年2月湖北省MODIS火点观测空间分布(a)以及与Himawari-8反演结果的散点图(b)Fig.4 The spatial distribution of fire events in Feb 2016 in Hubei province(a)and scatter plots of our FRP simulations and Himawari-8 FRP(b)

3.2 与地面环境监测站观测数据的对比结果

由于地面环境监测站观测的NO2浓度易受到来自交通源、工业源等非农业源排放的影响,为尽可能地减少除农业源外其他排放源的影响,提高地面站观测数据的抗鲁棒性,本研究在湖北省内选取了3个远离城市中心和交通干道的环境监测站点,即分别位于武汉市郊的农田(30.30 N,113.85 E)(简称武汉站)、十堰市郊的林区(32.39 N,111.04 E)(简称十堰站)、以及襄阳市郊的农田(32.09 N,112.21 E)(简称襄阳站),并将地面环境监测站的日均NO2浓度观测均值处理为月均NO2浓度,其中武汉站的时间跨度为2014年1月至2016年12月,襄阳站和十堰市站的时间跨度均为2015年1月至2016年12月。图5显示的是2014–2016年本研究估算的NOx排放量与3个地面站点观测的NO2浓度的比较。值得注意的是,本研究计算结果为NOx排放量,地面站观测数据为NO2浓度,而大气中NOx和NO2的比例较难准确测定。在大气导则中规定,NOx与NO2比例为4:3,然而实际情况则更为复杂,因此这里仅通过两者的相对变化评估本研究计算结果的准确性。由图5可知,本研究计算结果与武汉站和襄阳站的NO2实测浓度相关性较好,时间序列上变化趋势一致,相关系数分别为0.56和0.65。但十堰站NO2观测数据与本研究估算结果的相关性较弱,误差主要分布在2015年和2016年的1、2、3月,站点监测NO2浓度和本研究估算的NOx排放量在一定程度上呈现负相关。主要原因是由于该站点位于林区,冠层高度高,可能导致大量火点被遗漏。此外,NO、N2O等作为林地土壤中菌群与落叶等生物质产生反硝化作用的中间产物会向大气中释放,NO快速氧化后转化成为NO2,会对大气中NO2的浓度产生影响。综上,由于环境监测站实际观测的NO2浓度不仅与地表生物质燃烧排放有关,同时受到气象、大气物理化学过程以及其他排放源的影响,致使十堰市的监测结果与本研究结果差异较大。

图5 本研究结果和站点数据的对比。(a)武汉;(b)襄阳;(c)十堰Fig.5 Comparison of our monthly NOx emissions with field observations.(a)Wuhan;(b)Xiangyang;(c)Shiyan

3.3 与前人研究及已有数据库的对比结果

全球火灾排放数据库(GFED)主要提供全球的火灾排放信息[13]。本实验采用的数据源与GFED 4.1s产品均对亚像元级别的小火点进行了考量,因此两者的对比实验更具有客观性。由于农业源生物质燃烧排放数据在GFED数据库中仅提供1997年至2016年的数据,因此本研究主要对比2010–2016年间本研究估算结果与GFED数据库的结果。如图6所示,两者NOx排放量随时间的变化趋势整体一致,在2012–2014年间,农业源露天生物质燃烧NOx年排放量呈上升趋势,与GFED结果一致,本研究估算的农业源露天生物质燃烧NOx年排放量在2014年达到极大值,自2015年起排放量呈下降趋势,表明2015年我国颁布的秸秆燃烧禁令政策对空气污染控制起到了积极作用,而农业源露天生物质燃烧NOx年排放量上下波动,主要受种植结构、物候特征等多因素共同影响。而GFED数据库时序变化极大值与本研究实验结果存在细微的差异,造成这种差异的主要原因在于本研究基于像元地物类别先验概率的排放因子加权策略能更好地反应不同作物类别生物质燃烧排放。

根据已有文献结果,王书肖等[5]和彭立群等[35]的估算结果略低于本研究的估算结果,分别为47 Gg(2006年)和49 Gg(2009年),他们主要通过问卷调查的方式给定全国主产作物的秸秆产出和露天燃烧比率,以确定全国范围内各类秸秆经由露天燃烧消耗的质量,再进一步估算污染物的排放量。然而调查问卷难以准确表征整个区域的秸秆燃烧情况,并且未将不同区域不同秸秆类型差异对排放因子的贡献纳入考量。张晓荟[36]在问卷调查的基础上从地域和燃烧类型角度对各省的秸秆燃烧情况进行了分类,进一步细化了农业生物质燃烧污染物排放估算的过程,但仍无法避免调查问卷与实际燃烧之间的误差,其2017年的排放量为51 Gg。此外,为了表明本研究顾及火点像元类别先验概率的排放因子加权策略相比单一排放因子能更加精细化排放清单,本研究对比了采用两种不同排放因子设置策略的对比实验,如图6所示,其中虚线为单一排放因子,设置为3.11 g·kg-1,实线为混合加权排放因子。实验结果表明,采用单一排放因子计算结果远高于本研究计算结果、GFED数据库以及已有文献的结果。

图6 与GFED数据库、之前的研究以及单一排放因子的对比Fig.6 Comparison of this study with the GFED,previous studies and result of single emission factor

3.4 与OMI卫星数据的对比结果

为了探究农业源生物质露天燃烧排放对实际对流层中NOx空间分布的影响,将本研究排放清单的估算结果作为输入,利用WRF-Chem模型(Weather research and forecasting model coupled with chemistry)[37]在充分考虑垂直方向扩散和物种间物理化学作用的影响下,模拟了我国农业源生物质燃烧NOx排放空间分布。本研究选取农业源排放较严重的2015年10月和非农业源排放较严重的2016年2月,以比较农业源和非农业源排放对实际对流层NOx浓度的贡献。

对比发现2015年10月OMI NO2VCDs以及WRF-Chem NOx模拟结果相较2016年2月均显示出明显差异。如图7所示,2015年10月OMI观测在我国华北平原NO2对流层浓度明显高于2016年2月。其中河北南部的邯郸、石家庄以及河南安阳与山西省交接地带尤为明显,此外陕西西安周边地区也呈现出高NO2浓度水平。而华北平原作为我国重要的粮食基地,作物秸秆产量大,尚没有完全普及现代化秸秆处理手段。因此秋季秸秆焚烧仍然是当地高NO2浓度的主要因素。相较2015年10月,2016年2月我国NO2浓度水平较低,此时农事活动相比秋季不活跃,然而在我国东北地区,尤其是黑龙江,辽宁由于纬度较高,冬季气温在全国处于较低的水平,农村地区仍然保有秸秆焚烧以及燃煤取暖的现象,造成冬季东北地区的NO2浓度水平较秋季上升。

图7 OMI NO2 VCDs的结果(a)–(b)和WRF-Chem的模拟结果(c)–(d)Fig.7 OMI tropospheric NO2 VCDs(a)–(b)and WRF-Chem simulations(c)–(d)

在加入中国多尺度排放清单(MEIC)后,WRF-Chem模拟结果显示,2015年10月NOx对流层浓度在华北和华中地区略高于2016年2月,此外,2016年2月新疆乌鲁木齐NOx浓度呈现高值,主要由于高纬度地区冬季供暖造成化石燃料以及秸秆焚烧需求增加,引起NOx浓度上升,加之乌鲁木齐地处准格尔盆地南缘,地势低,不利于污染物扩散。综上,WRF-Chem模拟NOx对流层柱浓度结果与OMI NO2对流层柱浓度观测结果空间分布整体一致,呈现出中东部高,西部低的趋势。差异主要体现在WRF-Chem模拟在华中地区的NOx浓度高于OMI的NO2浓度水平,这可能是由于Aura卫星观测时间的限制(当地时间13:30),OMI观测难以准确地捕捉到我国地面农业源生物质燃烧峰值排放,因此对NOx排放强度存在低估。OMI观测与WRF-Chem模拟结果间存在较强的相关关系,均方根误差分别为0.54和0.49(图8),冬季存在更多的人为排放为排放清单增加了不确定性。

图8 WRF-Chem的模拟结果和OMI观测数据的对比。(a)2015年10月;(b)2016年2月Fig.8 Scatter plot of WRF-Chem simulations and OMI observations.(a)October 2015;(b)February 2016

4 结论

本研究基于多源数据,结合一种顾及火点像元地物类别先验概率的排放因子加权策略,通过利用高斯函数模拟地表生物质燃烧FRP昼夜变化,建立了我国2003–2020年逐日1 km田块尺度生物质燃烧NOx排放清单,并综合利用GFED数据库、已有的文献、Himawari-8观测数据、站点观测以及OMI NO2VCDs数据等,对本研究建立的排放清单进行评估。结果表明:本研究基于FRP昼夜周期分布假设模拟的FRP与Himawari-8 FRP一致性较好,两者的相关性达0.65。从时序上看,本研究结果与GFED数据库、地面环境观测站资料也展现出较好的一致性,除十堰站外,其余相关系数均大于0.5。我国农业源生物质燃烧NOx排放呈现冬春高、夏秋低的规律,农业源生物质燃烧NOx排放的峰值主要在12月至次年2月,谷值出现在5–7月。纵观2003–2020年农业源NOx排放,2003–2014年缓慢上升,2014–2016年由于政策法规的限制排放有所减少,2016–2020年排放有所回弹,说明需持续发挥政策对生物质燃烧的管控作用。从空间分布上,采用MEIC数据库中的电力源、民用源、工业源、交通源等非农业源清单结合本研究估算的排放清单,通过WRF-Chem模拟了2015年10月与2016年2月NO2对流层柱浓度,并与OMI NO2VCDs数据进行了比对,结果整体一致性较好。但鉴于农业源生物质燃烧排放相较于人为源排放较低,代表性略有不足。综上,基于遥感FRP数据可获得较为可靠的生物质燃烧大气污染物排放清单及其精细的空间分布,这对理解农业生产管理、大气污染控制等都具有重要意义。