CN111505016A - 一种秸秆焚烧排放清单估算方法 - Google Patents

Abstract

本发明公开一种秸秆焚烧排放清单估算方法，解决现有方法监测准确性差的问题。所述方法包含：获取VIIRS载荷L1B级数据，对每个像元计算亮度温度；根据波段M13、M16的亮度温度值先提取潜在火点，然后在所述潜在火点基础上提取绝对火点，再剔除非秸秆火点，得到秸秆焚烧点；对所述秸秆焚烧点计算不同类型秸秆的单个像元焚烧量；计算不同网格尺度的秸秆焚烧量，建立秸秆焚烧六种污染物排放清单。本发明实现了基于VIIRS载荷的秸秆焚烧排放清单制作。

Description

一种秸秆焚烧排放清单估算方法

技术领域

本发明涉及卫星遥感领域，尤其涉及一种秸秆焚烧排放清单估算方法。

背景技术

秸秆是指小麦、水稻、玉米、薯类、油料、棉花、甘蔗等农作物在收货籽实之后剩余的部分，秸秆焚烧不仅会造成严重的空气污染，还引起能见度下降给道路交通安全带来隐患。人为巡查的办法无法解决秸秆焚烧问题，有人利用视频系统来监控秸秆焚烧，这种监测方式在小范围可以及时掌握重点污染源的动态，提高检查的针对性和时效性，但毕竟费用太高。现有利用卫星遥感数据获取秸秆焚烧火点及其排放信息，可以节省人力物力，但未形成秸秆焚烧排放清单，无法识别秸秆类型和排放通量。

发明内容

本发明提供一种秸秆焚烧排放清单估算方法，解决现有方法监测准确性差的问题。

为解决上述问题，本发明是这样实现的：

本发明实施例指出一种秸秆焚烧排放清单估算方法，包含以下步骤：获取VIIRS载荷L1B级数据，对每个像元计算亮度温度；根据波段M13、M16的亮度温度值先提取潜在火点，然后在所述潜在火点基础上提取绝对火点，再剔除非秸秆火点，得到秸秆焚烧点；对所述秸秆焚烧点计算不同类型秸秆的单个像元焚烧量：

G_i,pixel＝S_i,pixel×p_i×f_i

其中，i为秸秆类型序号，pixel为像元序号，G_i,pixle为秸秆类型序号i、像元序号pixel的单个像元焚烧量，S_i,pixle为第i种类型秸秆的像元面积，p_i为第i种类型秸秆的农作物单产，f_i为第i种类型秸秆的秸秆产生系数；计算不同网格尺度的秸秆焚烧量，建立秸秆焚烧排放清单，所述不同网格尺度的秸秆焚烧量为：

其中，m为网格尺度，j为污染物种序号、j＝1～6对应的污染物种分别为PM、SO₂、NO_x、BC、OC、CO，

为第j种污染物种、m×m网格尺度的秸秆焚烧量，pixel为像元序号，δ_pixel为第pixel像元的秸秆焚烧像元标识，G_i,pixle为所述秸秆类型序号i、像元序号pixel的单个像元焚烧量，γ_i,j为秸秆类型序号i、污染物种序号j的秸秆焚烧排放因子。

进一步地，所述提取潜在火点的方法为：对每个像元的所述M13、M16的亮度温度值进行判断，若T_M13>A₁且T_M13-T_M16>A₂，则该像元为潜在火点，其中T_M13、T_M16分别为每个像元的所述M13、M16的亮度温度值，A₁、A₂分别为第一、第二阈值。

进一步地，所述在潜在火点基础上提取绝对火点的方法为：对每一个所述潜在火点设立一个包含该潜在火点的像元窗口，计算所述像元窗口内M13、M16和M13-M16波段的各像元亮度温度平均值和平均绝对偏差；

对所述潜在火点进行判断，若至少满足以下两个条件之一，则该潜在火点为所述绝对火点：条件一：

且σ_M13＞A₃；条件二：ΔT_M13,16＞A₄且T_M13＞A₅；其中，所述S₁、S₂、S₃分别为第一、第二、第三系数，A₃、A₄、A₅分别为第三、第四、第五阈值。

优选地，所述剔除非秸秆火点包括剔除工业热源的高温热异常点和森林火点。

优选地，所述方法还包含：剔除云像元。

优选地，在所述获取VIIRS载荷L1B级数据，对每个像元计算亮度温度的步骤之后，所述方法还包含：剔除水体像元。

优选地，所述VIIRS载荷L1B级数据为NASA官网数据或VIIRS载荷直接接收的数据经过辐射校正和几何校正。

进一步地，所述云像元的识别方法为：对每个像元计算波段I₁、I₂、I₃的表观反射率，波段M16、波段I₅的亮度温度值，若满足以下至少一个条件，则该像元为云像元：条件一：

且T_M16＜A₇；条件二：

且

条件三：

且

其中，

分别为所述每个像元计算波段I₁、I₂、I₃的表观反射率，T_M16、

分别为波段M16、I₅的亮度温度值，A₆～A₁₁分别为第六～第十一阈值。

进一步地，所述水体像元的判断方法为：对每个像元计算波段I₁、I₂、I₃的表观反射率，若至少满足以下两个条件之一，则该像元是水体像元：条件一：

且

条件二：

且

其中，

分别为所述每个像元计算波段I₁、I₂、I₃的表观反射率。

本发明有益效果包括：本发明提出的秸秆焚烧排放清单估算方法，利用NPP卫星VIIRS载荷的数据，对地面秸秆焚烧情况进行监测，以及对秸秆焚烧点的污染排放通量实现估算，本发明能快速、准确、全面建立秸秆排放清单，基于卫星数据可实现对全国范围甚至更大范围的秸秆焚烧情况监测和排放量的估算。

附图说明

此处所说明的附图用来提供对本发明的进一步理解，构成本发明的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1为一种秸秆焚烧排放清单估算方法流程实施例；

图2为一种包含水体像元剔除的秸秆焚烧排放清单估算方法流程实施例。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明具体实施例及相应的附图对本发明技术方案进行清楚、完整地描述。显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

秸秆焚烧不仅会造成严重的空气污染，还引起能见度下降给道路交通安全带来隐患。人为巡查的办法无法解决秸秆焚烧问题，有人利用视频系统来监控秸秆焚烧，这种监测方式在小范围可以及时掌握重点关注范围的动态，提高检查的针对性和时效性，但毕竟费用太高。为此，有人利用无人机开展秸秆禁烧巡查，其巡查的范围较视频相比要扩大了很多，但毕竟我国秸秆焚烧的面积太大，无人机仍无法满足秸秆焚烧监测的需要。早期厉青等人(2009)利用卫星遥感技术监测了2007年全国秸秆焚烧状况，卫星遥感技术可以较好的监测全国的秸秆焚烧态势，而且研究表明秸秆焚烧火点数变化趋势与空气污染指数具有较好的一致性，在污染物扩散条件不利时将导致空气质量明显下降。张思等(2019)使用CMAQ等模拟进行空气质量模拟，结合卫星火点数据、采用排放因子法进行污染物年排放量估算，模拟解释了2014年5月7日四川盆地内发生的一次油菜秸秆焚烧引起的重污染事件。在秸秆焚烧的污染物排放估算方面，王书肖等人采用自下而上的地面调查方法进行秸秆焚烧排放量估算(王书肖等，2008；李建峰，2015)。针对这钟估算方法，彭立群等研究指出秸秆焚烧PM2.5排放量不确定性在-61％～99％之间(彭立群等，2016)。为此，徐敬等人(2018)采用卫星监测等火点燃手排放数据，利用区域化学传输模式WRF-Chem模拟分析了2017年5月秸秆焚烧现象对华北地区空气质量照成一定影响天数达20天。毛慧琴等人(2018)利用MODIS卫星遥感数据获取的秸秆焚烧火点与站点监测信息进行空间分析，即能大量节省人力物力，又有很好的时效性，并揭示其时空分布特性。但到目前为止，还缺少利用卫星遥感估算秸秆焚烧排放清单的方法，使得在区域空气质量预报中，往往由于缺乏秸秆焚烧排放清单的及时输入，造成空气质量预报结果出现巨大的误差。所以，本提案基于现有的卫星遥感监测秸秆焚烧技术，提出快速估算秸秆焚烧污染物排放清单技术。

本发明创新点如下：第一，目前还没有针对VIIRS载荷数据实现秸秆焚烧点监测的方法，本发明创新性地针对VIIRS载荷数据实现了秸秆焚烧点监测；第二，本发明创新性地对秸秆焚烧产生的PM、SO₂、NO_x、BC、OC、CO六种污染物种实现了排放通量估算，细化了秸秆焚烧排放清单，使得清单更准确全面。

以下结合附图，详细说明本发明各实施例提供的技术方案。

图1为一种秸秆焚烧排放清单估算方法流程实施例，利用卫星遥感数据对地面秸秆焚烧点进行监测，作为本发明实施例，一种秸秆焚烧排放清单估算方法，具体包含以下步骤：

步骤101，获取VIIRS载荷L1B级数据，对每个像元计算亮度温度。

在步骤101中，对每个像元可见光-近红外和短波红外波段计算表观反射率，对每个像元中红外-热红外波段计算亮度温度。

在步骤101中，利用美国NPP卫星VIIRS载荷数据进行秸秆焚烧监测，所述VIIRS载荷L1B级数据为NASA官网数据或VIIRS载荷直接接收的数据经过辐射校正和几何校正。

NASA(National Aeronautics and Space Administration，美国航空航天局)正在构建NPOESS Preparatory Project卫星系统(简称NPP)，由NASA、NOAA(NationalOceanic and Atmospheric Administration，美国国家海洋和大气管理局)和美国空军共同研发。2011年NASA发射了首颗NPP卫星，星上搭载了可见光红外成像辐射仪VIIRS(Visible infrared Imaging Radiometer)，是基于高分辨率辐射仪AVHRR(Advanced VeryHigh Resolution Radiometer)和中分辨率成像光谱仪MODIS(Moderate-resolutionImaging Spectroradiometer)两个载荷中优选波段设计的，是新一代极轨探测载荷。

在步骤101中，所述VIIRS载荷L1B级数据可以为NASA官网数据，从NASA的卫星数据网站直接下载VIIRS载荷L1B级数据，该数据已经进行了辐射校正和几何校正，VIIRS L1B下载地址：

https://ladsweb.modaps.eosdis.nasa.gov/search/

在步骤101中，所述VIIRS载荷L1B级数据为利用自己接收的VIIRS数据，需要对数据进行辐射校正和几何校正。需要说明的是，辐射校正和几何校正是卫星遥感数据定量化应用的两个基本的成熟的处理技术，辐射校正主要利用辐射校正模型：

L_b0,b1＝Gain_b0×DN_b0,b1+Offset_b0 (1)

其中，DN_b0,b1、L_b0,b1分别为第b0波段第b1像元的整型数值和经过校正的辐射值，Gain_b0、Offset_b0分别为第b0波段的第一、第二定标系数，从NASA的数据集中可以读取。

几何校正是指消除或改正遥感影像几何形变误差的过程。遥感影像的变形误差，大体分为两类：静态误差和动态误差，静态误差是在成像过程中，传感器相对与地球表面呈静止状态时所具有的各种形变误差；动态误差主要是在成像过程中由于地球旋转等因素造成的图像变形误差，而变形误差又可分为内部误差和外部误差两类。具体的处理算法已经在威斯康辛大学的卫星遥感数据定量处理软件包IPOPP(International PlanetaryObservation Processing Package)中实现，可以利用该软件进行处理。该软件的下载地址：

http://dbmeeting.sci.gsfc.nasa.gov/。

步骤102，根据波段M13、M16的亮度温度值先提取潜在火点，然后在所述潜在火点基础上提取绝对火点，再剔除非秸秆火点，得到秸秆焚烧点。

在步骤102中，所述提取潜在火点的方法为：对每个像元的所述M13、M16的亮度温度值进行判断，若T_M13>A₁且T_M13-T_M16>A₂，则该像元为潜在火点，其中T_M13、T_M16分别为每个像元的所述M13、M16的亮度温度值，A₁、A₂分别为第一、第二阈值。

需要说明的是，所述第一、第二阈值为设定数值，这里不做特别限定。

例如，根据中红外波段对高温异常探测的敏感性，利用VIIRS载荷波段M13：3.9735～4.128μm和波段M6：11.538～12.488μm的亮度温度，以及中红外和热红外波段的亮温差来确定所述潜在火点，确定方式为：T_M13>290K且T_M13-T_M6>10K，即所述第一阈值A1＝290K，所述第二阈值A2＝10K。

在步骤102中，所述在潜在火点基础上提取绝对火点的方法为：对每一个所述潜在火点设立一个包含该潜在火点的像元窗口，如窗口尺度为m×m，计算所述像元窗口内M13、M16和M13-M16波段的各像元亮度温度平均值和平均绝对偏差：

其中，m为所述像元窗口的网格尺度，i₀为包含在所述像元窗口内的像元序号，

分别为所述M13、M16、M13-M16波段的亮度温度平均值，σ_M13、σ_M16、σ_M13，16分别为所述M13、M16、M13-M16波段的平均绝对偏差，

ΔT_i0,M13,16分别为所述M13、M16、M13-M16波段、像元序号为i₀的亮度温度值。

对所述潜在火点进行判断，若至少满足以下两个条件之一，则该潜在火点为所述绝对火点：

条件一：

且σ_M13＞A₃；

条件二：ΔT_M13,16＞A₄且T_M13＞A₅；

其中，所述S₁、S₂、S₃分别为第一、第二、第三系数，A₃、A₄、A₅分别为第三、第四、第五阈值。

需要说明的是，m₁、所述第一系数、第二系数、第三系数、第三阈值、第四阈值、第五阈值均为设定数值，这里不做特别限定。

例如，针对潜在火点探测结果，设立一个“21×21”像元窗口来比较潜在火点和窗口背景温度的状况，即所述像元窗口的尺度大小m₁为21。先计算VIIRS载荷波段M13、M16和两个波段差在窗口内各像元亮度温度的平均值和平均绝对偏差，当像元满足条件一：

且σ_M13＞5K，或满足条件二：ΔT_M13,16＞25K且T_M13＞330K，则像元为所述绝对火点，即所述第一系数S₁为3、第二系数S₂为1、第三系数S₃为3.5，所述第三阈值A₃为5K、第四阈值A₄为25K、第五阈值A₅为330K。

在步骤102中，所述剔除非秸秆火点包括剔除工业热源的高温热异常点和森林火点。

具体地，在探测出的火点中，实际上有一类是工业热源的高温异常点，工业热源也会被检测到。利用2018年监测到的高温异常数据，把位于1km²范围之内的高温异常点剔除，因为工业热源的位置不变，但考虑其几何定位有半个像元的误差，所以同一个工业热源会出现热点团簇现象，将这些点建立一个工业热源数据库，把今后监测到的高温异常点且位于该范围的，归属于工业热源。

需要说明的是，除利用2018年监测到的高温异常数据外、也可利用其它高温异常数据，距离范围可以是1km²、也可以是其他数值，这里均不作特别限定。

另外，由于真实火点中有森林火点现象，所以需要结合清华大学土地利用与土地覆盖数据(http://data.ess.tsinghua.edu.cn)，把位于农田范围内的火点提取出来，作为秸秆焚烧火点信息。

需要说明的是，本发明结合清华大学土地利用与土地覆盖数据，也可结合其他土地覆盖和土地利用数据来去除森林火点，这里不做特别限定。

步骤103，对所述秸秆焚烧点计算不同类型秸秆的单个像元焚烧量：

G_i,pixel＝S_i,pixel×p_i×f_i (9)

其中，i为秸秆类型序号，pixel为像元序号，G_i,pixle为秸秆类型序号i、像元序号pixel的单个像元焚烧量，S_i,pixle为第i种类型秸秆的像元面积(ha)，p_i为第i种类型秸秆的农作物单产，f_i为第i种类型秸秆的秸秆产生系数。

在步骤103中，针对秸秆焚烧，农民不会收集秸秆一块堆烧，而是点着秸秆让其在地里焚烧蔓延，这样导致卫星探测到的像元有秸秆焚烧点后，该像元秸秆按烧尽处理。

设三种主要农作物单产为p_i，单位是“公斤/公顷(kg/ha)”，在本发明中i＝1,2,3分别代表小麦、水稻和玉米。这个数据可以从国家相关部门公布的数据中获得，如根据新浪财经报道，我国2019年全国粮食单位面积产量5720公斤/公顷(新浪财经，2019年12月10日)；根据全国第一次污染源普查农业源普查数据，和产量相比，小麦、水稻和玉米的秸秆产生系数分别是小麦f₁＝1.2、水稻f₂＝1.06、玉米f₃＝1.34，秸秆产生系数为无量纲。需要说明的是，农作物单产可以通过2019年数据得到、也可通过其他数据得到，这里不做特别限定；所述秸秆产生系数，可以通过全国第一次污染源普查农业源普查数据得到本实施例数值、也可通过其他方式得到其他数值，这里不做特别限定。

在步骤103中，计算卫星影像每个像元的秸秆产量，由于美国NASA卫星VIIRS载荷的中红外和热红外波段的空间分辨率正好是1km，所以每个像元的面积正好是：1km×1km＝100ha，这样每个像元的秸秆产量模型是“像元面积×单位面积粮食单产×不同粮食秸秆系数”。

步骤104，计算不同网格尺度的秸秆焚烧量，建立秸秆焚烧排放清单，所述不同网格尺度的秸秆焚烧量为：

其中，m为网格尺度，j为污染物种序号、j＝1～6对应的污染物种分别为PM、SO₂、NO_x、BC、OC、CO，

为秸秆类型序号i、第j种污染物种、m×m网格尺度的秸秆焚烧量，pixel为像元序号，δ_pixel为第pixel像元的秸秆焚烧像元标识，G_i,pixle为所述秸秆类型序号i、像元序号pixel的单个像元焚烧量，γ_i,j为秸秆类型序号i、污染物种序号j的秸秆焚烧排放因子。

需要说明的是，若像元是秸秆焚烧像元，则pixel＝1；否则，若像元是非秸秆焚烧像元，则pixel＝0。

还需说明的是，m在步骤104中的与步骤102中相同，均表示像元窗口的网格尺度。

在步骤104中，计算所述秸秆焚烧排放因子的方法为现有技术，通过实验室数据获得，可以采用本发明实施例中的计算方法和数值，也可采用其他计算方法和数值，这里不做特别限定。

在本发明实施例中，秸秆焚烧排放污染物的量和秸秆的有机质含量、含水率、燃烧温度和环境风速、温度等自然环境状况有关，根据自全国各个有代表性的秸秆进行实验室燃烧实验，包括麦秸秆、稻秸秆、玉米杆、棉花杆，实验测量了PM、SO₂、NO_x、BC、OC和CO的排放因子(曹国良等，2007)，具体不同类型秸秆的排放因子为：

i＝1代表小麦秸秆焚烧排放因子(g/kg)，j＝1代表污染物种PM，则小麦秸秆PM的焚烧排放因子γ_1,1＝9.64，同理，小麦秸秆SO₂的焚烧排放因子γ_1,2＝0.049，小麦秸秆NO_x的焚烧排放因子γ_1,3＝2.59，小麦秸秆BC的焚烧排放因子γ_1,4＝0.52，小麦秸秆OC的焚烧排放因子γ_1,5＝3.83，小麦秸秆CO的焚烧排放因子γ_1,6＝65.5。

同样地，i＝2代表水稻秸秆焚烧排放因子，水稻秸秆PM的焚烧排放因子γ_2,1＝6.04，同理，水稻秸秆SO₂的焚烧排放因子γ_2,2＝0.147，水稻秸秆NO_x的焚烧排放因子γ_2,3＝3.52，水稻秸秆BC的焚烧排放因子γ_2,4＝0.52，水稻秸秆OC的焚烧排放因子γ_2,5＝1.96，水稻秸秆CO的焚烧排放因子γ_2,6＝72.4。

同样地，i＝3代表玉米秸秆焚烧排放因子，玉米秸秆PM的焚烧排放因子γ_3,1＝5.26，同理，玉米秸秆SO₂的焚烧排放因子γ_3,2＝0.026，玉米秸秆NO_x的焚烧排放因子γ_3,3＝3.36，玉米秸秆BC的焚烧排放因子γ_3,4＝0.78，玉米秸秆OC的焚烧排放因子γ_3,5＝2.21，玉米秸秆CO的焚烧排放因子γ_3,6＝70.2。

在步骤104中，由于区域大气化学模式的网格尺度根据需求分为不同层次，其要求排放清单数据也需要按相匹配的网格尺度的文档存储。所以，针对不同网格尺度需要分别建立清单。

一般地，所述网格尺度为8km或64km，需要说明的是，所述网格尺度也可以其他数值，这里不做特别限定。

在步骤104中，单个像元污染物排放通量估算为，针对窗口每个像元的情况，其排放量估算需要先区分秸秆焚烧像元和非秸秆焚烧像元，其估算模型为：

其中，

为单个像元、第i类秸秆、第j种污染物种的秸秆焚烧量。

在步骤104中，针对窗口尺度m，则m×m个像元的秸秆焚烧量如公式10，分别计算了小麦秸秆的PM、SO₂、NO_x、BC、OC、CO六种污染物种的秸秆焚烧量，即污染物排放通量；计算了水稻秸秆的PM、SO₂、NO_x、BC、OC、CO六种污染物种的秸秆焚烧量；计算了玉米秸秆的PM、SO₂、NO_x、BC、OC、CO六种污染物种的秸秆焚烧量。

需要说明的是，本发明方法适用的秸秆类型不限于小麦、玉米、水稻，还适用于其他秸秆类型。

本发明实施例基于VIIRS载荷数据，根据秸秆类型、排放污染物种类型建立了排放清单并计算了排放通量，细化了秸秆类型、排放污染物，使得秸秆排放清单更准确、全面，更具有实用价值。

图2为一种包含水体像元剔除的秸秆焚烧排放清单估算方法流程实施例，作为本发明实施例，一种秸秆焚烧排放清单估算方法，具体包含以下步骤：

步骤201，获取VIIRS载荷L1B级数据，对每个像元计算亮度温度。

步骤201与步骤101相同，这里不做具体说明。

步骤202，剔除水体像元和云像元。

在步骤202中，VIIRS L1B数据在0.3～14μm波段范围内有22个，其中可见光、近红外9个(0.4～0.9μm)；短、中波红外8个(1～4μm)；热红外4个(8～12μm)；秸秆焚烧的监测算法需要应用可见光-近红外波段和短波红外波段，以及中红外和热红波段，为此需要相关波段进行表观发射率与亮度温度的数据处理。

对可见光-近红外和短波红外波段的数据进行表观反射率计算，具体利用下面公式计算：

其中，ρ为大气顶表观反射率(无量纲)，E_sun为大气顶层太阳辐照度，单位W/(m²srμm)，θ为太阳天顶角，L为大气顶入瞳波段辐亮度，单位W/(m²sr μm)，D为日地之间距离(天文单位)，sr为单位立体角。

针对中红外和热红外波段，需要计算将每个波段的辐射值转换为亮度温度值，具体基于下面模式计算：

其中，这里λ是某波段的中心波长值，单位μm，T是亮度温度，单位是绝对温度K，L₀是像元某波段的辐射亮度值，单位W/(m²sr μm)；c是真空中的光速；K是玻耳兹曼常数1.38×10^-23J/K，h是普朗克常数，6.626×10^-34Js。

在步骤202中，云层能阻挡地表秸秆焚烧热辐射信息的传输，所以要首先开展对云覆盖像元的识别，依据VIIRS载荷的波段设置，进一步地，所述云像元的识别方法为：

对每个像元计算波段I₁、I₂、I₃的表观反射率，波段M16、波段I₅的亮度温度值，若满足以下至少一个条件，则该像元为云像元：

条件一：

且T_M16＜A₇；

条件二：

且

条件三：

且

其中，

分别为所述每个像元计算波段I₁、I₂、I₃的表观反射率，T_M16、

分别为波段M16、I₅的亮度温度值，A₆～A₁₁分别为第六～第十一阈值。

需要说明的是，A₆为所述第六阈值、A₇为所述第七阈值、A₈为所述第八阈值、A₉为所述第九阈值、A₁₀为所述第十阈值、A₁₁为所述第十一阈值，均为可设定数值，对具体数值这里不做特别限定。

例如，探测条件一：

且T_M16＜265K；探测条件二：

且T_I5＜265K；探测条件三：

且T_I5＜285K，其中，

分别为所述每个像元计算波段I₁、I₂、I₃的表观反射率，T_M16、T_I5分别为波段M16、I5的亮度温度值，单位K，所述第六阈值为0.9，所述第七阈值为265K，所述第八阈值为0.6，所述第九阈值为265K，所述第十阈值为0.4，所述第十一阈值为285K。

需要说明的是，VIIRS载荷中有两种波段序号，一种是以“I”为标志的，如所述的I1、I2、I3等，分别有波段I₁：0.600～0.680μm，波段I₂：0.846～0.885μm，波段I₃：1.580～1.640μm，波段I₅：10.500～12.400μm；一种序号是以“M”为标志的，如所述的M13、M16等，有M13：3.9735～4.128μm和波段M16：11.538～12.488μm。

在步骤202中，针对水体表面，由于卫星载荷的观测与太阳观测几何形成镜面反射时，会产生镜面反射现象而导致假火点的产生，所以在秸秆焚烧监测算法中，我们将卫星影像中的水体覆盖像元全部去除，具体算法主要是依据水体的近红外反射率非常小的特点来开展识别。具体算法如下：

对每个像元计算波段I₁、I₂、I₃的表观反射率，若至少满足以下两个条件之一，则该像元是水体像元：

条件一：

且

条件二：

且

其中，

分别为所述每个像元计算波段I₁、I₂、I₃的表观反射率，A₁₂为第十二阈值。

需要说明的是，所述第十二阈值为设定数值，这里不做具体限定。

例如，条件一：

且

条件二：

且

其中，

分别为所述每个像元计算波段I₁、I₂、I₃的表观反射率，所述第十二阈值为0.1。

步骤203，根据波段M13、M16的亮度温度值先提取潜在火点，然后在所述潜在火点基础上提取绝对火点，再剔除非秸秆火点，得到秸秆焚烧点。

在步骤203中，提取潜在火点的数据为步骤202中经过剔除水体像元的数据，得到所述秸秆焚烧点的方法与步骤102相同。

步骤204，对所述秸秆焚烧点计算不同类型秸秆的单个像元焚烧量。

步骤204计算所述单个像元焚烧量的方法与步骤103相同。

步骤205，计算不同网格尺度的秸秆焚烧量，在所述秸秆焚烧排放清单中加入观测时间，建立秸秆焚烧排放清单。

需要说明的是，本发明实施例建立了六种污染物排放清单，六种污染物分别为PM、SO₂、NO_x、BC、OC、CO。

在步骤205中，计算不同网格尺度的秸秆焚烧量的方法与步骤104相同，这里不重复论述。

需要说明的是，在所述秸秆焚烧排放清单中加入观测时间可以说明是日监测还是月尺度等不同时间尺度的观测。

在步骤205中，对不同网格尺度的秸秆焚烧量数据存储好，并标出观测时间，也可根据一定时间尺度对每天排放量的累计计算以实现秸秆燃烧排放通量文档制作，以满足不同时间和空间分辨率模拟需求。

本发明实施例提供的秸秆焚烧排放清单估算方法，剔除水体像元影响，使焚烧点清单更准确，另外，针对不同网格尺度的秸秆焚烧量标记了时间标签，使数据更具有可塑性，工程应用型更强。

需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、商品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、商品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、商品或者设备中还存在另外的相同要素。

以上所述仅为本发明的实施例而已，并不用于限制本发明。对于本领域技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (9)

1.一种秸秆焚烧排放清单估算方法，其特征在于，包含以下步骤：

获取VIIRS载荷L1B级数据，对每个像元计算亮度温度；

根据波段M13、M16的亮度温度值先提取潜在火点，然后在所述潜在火点基础上提取绝对火点，再剔除非秸秆火点，得到秸秆焚烧点；

对所述秸秆焚烧点计算不同类型秸秆的单个像元焚烧量：

G_i,pixel＝S_i,pixel×p_i×f_i

其中，i为秸秆类型序号，pixel为像元序号，G_i,pixle为秸秆类型序号i、像元序号pixel的单个像元焚烧量，S_i,pixle为第i种类型秸秆的像元面积，p_i为第i种类型秸秆的农作物单产，f_i为第i种类型秸秆的秸秆产生系数；

计算不同网格尺度的秸秆焚烧量，建立秸秆焚烧排放清单，所述不同网格尺度的秸秆焚烧量为：

其中，m为网格尺度，j为污染物种序号、j＝1～6对应的污染物种分别为PM、SO₂、NO_x、BC、OC、CO，

为秸秆类型序号i、第j种污染物种、m×m网格尺度的秸秆焚烧量，pixel为像元序号，δ_pixel为第pixel像元的秸秆焚烧像元标识，G_i,pixle为所述秸秆类型序号i、像元序号pixel的单个像元焚烧量，γ_i,j为秸秆类型序号i、污染物种序号j的秸秆焚烧排放因子。

2.如权利要求1所述的秸秆焚烧排放清单估算方法，其特征在于，所述提取潜在火点的方法为：

对每个像元的所述M13、M16的亮度温度值进行判断，若T_M13>A₁且T_M13-T_M16>A₂，则该像元为潜在火点，其中T_M13、T_M16分别为每个像元的所述M13、M16的亮度温度值，A₁、A₂分别为第一、第二阈值。

3.如权利要求1所述的秸秆焚烧排放清单估算方法，其特征在于，所述在潜在火点基础上提取绝对火点的方法为：

对每一个所述潜在火点设立一个包含该潜在火点的像元窗口，计算所述像元窗口内M13、M16和M13-M16波段的各像元亮度温度平均值和平均绝对偏差：

其中，m为所述像元窗口的网格尺度，i₀为包含在所述像元窗口内的像元序号，

分别为所述M13、M16、M13-M16波段的亮度温度平均值，σ_M13、σ_M16、σ_M13，16分别为所述M13、M16、M13-M16波段的平均绝对偏差，

分别为所述M13、M16、M13-M16波段、像元序号为i₀的亮度温度值；

对所述潜在火点进行判断，若至少满足以下两个条件之一，则该潜在火点为所述绝对火点：

条件一：

且σ_M13＞A₃；

条件二：ΔT_M13,16＞A₄且T_M13＞A₅；

其中，所述S₁、S₂、S₃分别为第一、第二、第三系数，A₃、A₄、A₅分别为第三、第四、第五阈值。

4.如权利要求1所述的秸秆焚烧排放清单估算方法，其特征在于，所述剔除非秸秆火点包括剔除工业热源的高温热异常点和森林火点。

5.如权利要求1所述的秸秆焚烧排放清单估算方法，其特征在于，所述方法还包含：剔除云像元。

6.如权利要求1所述的秸秆焚烧排放清单估算方法，其特征在于，在所述获取VIIRS载荷L1B级数据，对每个像元计算亮度温度的步骤之后，所述方法还包含：剔除水体像元。

7.如权利要求1所述的秸秆焚烧排放清单估算方法，其特征在于，所述VIIRS载荷L1B级数据为NASA官网数据或VIIRS载荷直接接收的数据经过辐射校正和几何校正。

8.如权利要求5所述的秸秆焚烧排放清单估算方法，其特征在于，所述云像元的识别方法为：

对每个像元计算波段I₁、I₂、I₃的表观反射率，波段M16、波段I₅的亮度温度值，若满足以下至少一个条件，则该像元为云像元：

条件一：

且T_M16＜A₇；

条件二：

且

条件三：

且

其中，

分别为所述每个像元计算波段I₁、I₂、I₃的表观反射率，T_M16、

分别为波段M16、I₅的亮度温度值，A₆～A₁₁分别为第六～第十一阈值。

9.如权利要求6所述的秸秆焚烧排放清单估算方法，其特征在于，所述水体像元的判断方法为：

对每个像元计算波段I₁、I₂、I₃的表观反射率，若至少满足以下两个条件之一，则该像元是水体像元：

条件一：

且

条件二：

且

其中，

分别为所述每个像元计算波段I₁、I₂、I₃的表观反射率，A₁₂为第十二阈值。

Images