ENVI下的秸秆焚烧实时遥感监测方案及实现
分享
-
2016-06-24
每年夏收与秋收后，部分城市始终笼罩着浓浓的烟雾。这些烟雾不是来自重工业污染，也不是化工产品爆炸，而是来自秸秆焚烧。秸秆焚烧已经成为影响这些城市空气质量的重要因素，不仅如此，持续大雾还会使得重点城市的交通运输效率降低，甚至引发交通事故。秸秆焚烧严重影响大气环境质量，导致空气中总悬浮颗粒物数量明显升高，而焚烧产生的浓烟中含有大量的CO、CO<sub>2</sub>等气体，这样的气体刺激呼吸道，对人体健康产生不良影响。
图1：2012年6月11号的武汉大雾
鉴于秸秆焚烧带来的种种危害，秸秆焚烧的监测已经引起了各级人民政府的高度关注，利用实时监测结果并依照相关法律法规对其即时制止是杜绝秸秆焚烧的首要手段。而传统的监测手段（如逐点人工排查）具有效率低、覆盖率低、成本高等缺点。
卫星遥感手段以其时效性、覆盖面广、分辨率高等优势使得快速大面积监测焚烧情况成为可能。MODIS是先进的多光谱遥感传感器，具有36个观测通道，覆盖了当前主要遥感卫星的主要观测数据。利用MODIS的陆地产品提取秸秆焚烧火点的处理流程如下：
图2：技术流程图
流程涉及到影像的光谱处理和几何处理。ENVI/IDL遥感影像处理软件能对MODIS数据进行快速精确的光谱和几何处理，给更快速、更准确的监测秸秆焚烧现状提供了技术保障。其中大气校正模块FLAASH是采用精度最高大气校正算法，是基于光谱信息提取相关专题产品的必要前提；其面向对象信息提取工具ENVI EX能有效的对地表进行分类，提取土地覆盖信息；其强大的波段运算工具是提取各种专题数据的有效手段。
由于MODIS的 MOD14热异常数据可供直接获取使用，能够探测比气象卫星更小更多的火点（面积50平方米），是监测秸秆焚烧理想的数据源。下面从数据源获取、火点提取、农用地提取、结果输出等几个方面来介绍ENVI在秸秆焚烧监测中的应用方案。
1、数据获取
MODIS数据下载地址：http://ladsweb.nascom.nasa.gov/data/search.html。
(一) 首先打开网页，选中需要的数据源类型。其中，日期类型为：月/日/年 时:分:秒，其中，网页中显示的时间为UTC时间。
(二) 在“spatial selection” 选项中选择“latitude/longtitude”，按监测区域的经纬度选择影像范围。点击“search”查到需要的数据，勾选需要的数据，点击“order files now”， 输入你接收信息的邮箱，点“order”开始订购该数据。(如果要搜索多天数据，可以选“add files to shopping cart”继续搜索其他日期的数据。
(三) 所订购数据的存放位置信息： 点击“Data->Track Orders ”可以查看所有已订购的数据的状态。如果“state”显示“avalable”即可开始下载。
本实验数据覆盖范围：
图3：数据覆盖范围示意图
2、火点提取
图4：Modis数据
说明：在这个步骤中，一般情况下，如果不是采用MOD14数据，而是由MOD02数据以及MOD03数据进行热异常信息提取，则需要用到ENVI的FLAASH大气校正模块进行大气校正。在这里，我们简化地采用MOD14数据。
第一步：MOD14数据的地理定位
MOD14热异常数据地理定位必须要用到MOD03数据。
1、 “MOD14.A2012163.0305.005.2012163094359.hdf”显示如下：
图5：MOD14热异常数据集
2、 选择MOD03的经纬度数据打开：
图6：MOD03数据集信息
3、利用经纬度构建GLT地理查找表
图7：GLT文件构建
“X Geometry Band”即为“Longitude”， “Y Geometry Band”即为“Latitude”：选择输出坐标系,设置GLT参数，创建完成后，在Available Band List出现。
图8：GLT文件构建
4、利用GLT地理查找表对MOD14数据进行地理定位。
利用上一步生成的GLT地理查找表对MOD14进行地理定位，坐标系为WGS-84，结果如下：
图9：地理定位后的MOD14热异常数据
第二步：火点提取
对MOD14地理定位后的波段进行波段信息统计，统计信息如下：
图10：MOD14统计结果
0-9分别代表的含义如下表所示：
图11：MOD14数据像元值分类
可见，像素值为7，8，9的地方为火点。利用波段运算窗口进行火点提取：
BandMath表达式：(b1 ge 7)*0+(b1 lt 7)*1
注：就是讲7、8、9像元值变成1，其余为0，对MOD14数据进行二值化处理。
图12：提取火点信息的BandMath表达式
在表达式窗口中输入二值化表达式对图像进行二值化，火点区赋值为0，非火点区赋值为1. 得到的结果如下图所示：
图13：MOD14二值化结果
利用感兴趣区域将火点影像转化为矢量，结果（红色三角形点）如下：
图14：疑似火点矢量化结果
3、农用地提取
说明：在此步骤中，一般情况下，需要使用实时 MOD02数据进行农用地提取，此时需要用到ENVI EX面向对象信息提取工具对土地覆盖进行分类提取。在这里我们简华地采用现势性不强（2009年）的MCD12Q1数据作为替代。
火点提取出来后，接下来需要判定火点是否在农用地之上，在此，我们认为在农用地之上的火点为秸秆焚烧造成的火点。农用地的提取方法有多种，在此我们重点介绍由MODIS三级数据土地覆盖类型产品提取农用地的方法，其余方法则不一一列举。
MODIS三级数据土地覆盖类型产品（Land Cover data）是根据一年的Terra和Aqua观测所得的数据经过处理，描述土地覆盖的类型。该土地覆盖数据集中包含了17个主要土地覆盖类型，根据国际地圈生物圈计划（IGBP），其中包括11个自然植被类型，3个土地开发和镶嵌的地类和3个非草木土地类型定义类。
MODIS Terra + Aqua三级土地覆盖类型年度全球500米产品MCD12Q1采用五种不同的土地覆盖分类方案，信息提取主要技术是监督决策树分类。下面是该数据中包含的五个数据集，五个分类方案如下：
土地覆盖分类1：IGBP的全球植被分类方案
土地覆盖分类2：美国马里兰大学（UMD格式）方案
土地覆盖分类3：基于MODIS叶面积指数/光合有效辐射方案
土地覆盖分类4：基于MODIS衍生净初级生产力（NPP）方案
土地覆盖分类5：植物功能型（肺功能）方案
MCD12Q1是HDF格式的，包含16个数据集
图15：MCD12Q1数据集说明
第一步：数据读取
我们利用MCD12的第一类分类数据进行农用地信息提取。
图16：MCD12Q1数据集信息
（1）Land_cover_type数据说明
打开Land_cover_type1到Land_cover_type5，是五类分类方案的分类图，图像上的DN值各个值都代表某一地物类型，具体意义如下表：
图17：MCD12Q1数据分类说明
我们感兴趣的是第一种分类图像中DN值为12的像元。
第二步：数据地理定位
MODIS三级数据产品土地覆盖数据是已经经过投影了的，其投影类型为正弦曲线投影，但是很多遥感软件在打开HDF数据的时候，会有投影丢失现象，ENVI也是，可以通过hdf属性查看工具，查看ENVI打开hdf时未识别的投影信息，以及投影参数，利用该参数在ENVI中自定义正弦曲线投影，为MCD12Q1数据添加投影。
（1）用hdf属性查看工具查看该hdf数据的地理定位信息，可得图像左上角、右下角坐标值，投影类型，投影参数，地球半径，起始位置。
图18：MCD12Q1数据集的地理定位信息
（2）根据投影信息，在ENVI中自定义投影：
图19：自定义正弦曲线投影坐标
并修改其头文件信息：
图20：编辑MapInfo坐标信息MCD12Q1数据
已经具有投影信息的数据：
图21：定义Mapinfo坐标信息后的MCD12Q1数据
第三步：投影变换与镶嵌
（1） 针对每一张MODIS土地覆盖的影像，经过地理定位后，将其坐标转到WGS-84坐标系下，用投影转换工具对相应的影像进行投影转换，变换到WGS84坐标系下：
图22：MCD12Q1数据投影转换
（2）影像镶嵌
所有的土地覆盖影像进行投影转换之后，将此数十幅影像进行拼接。用基于地理投影的镶嵌工具进行影像镶嵌，如下图所示：
图23：MCD12Q1数据镶嵌面板
将所有的背景值设为0，避免互相压盖。镶嵌结果如下：
图24：MCD12Q1数据镶嵌结果
同样，利用火点提取一节中的方法，依照MCDQ1的分类说明，把DN=12的像素（农耕用地）提取出来。结果如下图所示（绿色区域）
图25：从MCD12Q1数据集中提取农耕用地
4、结果输出
在Image窗口中将结果连接至ArcMap，第二个环节的火点与第三环节的土地覆盖叠加显示如下：
图26：火点提取结果
该结果可以在ArcGIS中制图制图输出或发布服务。
5、IDL定制开发
ENVI是基于IDL语言开发的，也提供了丰富的开发接口供用户调用，上述流程在ENVI中都有相应的API，利用IDL语言将各个流程串接起来，可以形成专门的工具嵌入ENVI主界面，或作为一个独立的IDL程序提供给软件使用者，亦可嵌入其他语言做系统集成开发，可以降低对遥感专业知识的要求，直接从原始数据到最终结果。当然，开发需要专业的开发人员参与。
下面是基于ENVI/IDL开发的基于HJ的火点检测的程序（2009中南大学，彭光雄），火点检测采用如下算法（结果与MOD14相当）：
程序嵌入在ENVI主菜单中：
此时，只需要点击菜单，选择影像：
输入少数参数：
就可以得到火点分布情况：
该程序可以作为ENVI的一个插件，也可以作为一个单独的应用程序，或是与C# /Java做混合开发嵌入业务系统中。
文章来源：http://blog.sina.com.cn/s/blog_764b1e9d01010sl9.html
图1：2012年6月11号的武汉大雾
鉴于秸秆焚烧带来的种种危害，秸秆焚烧的监测已经引起了各级人民政府的高度关注，利用实时监测结果并依照相关法律法规对其即时制止是杜绝秸秆焚烧的首要手段。而传统的监测手段（如逐点人工排查）具有效率低、覆盖率低、成本高等缺点。
卫星遥感手段以其时效性、覆盖面广、分辨率高等优势使得快速大面积监测焚烧情况成为可能。MODIS是先进的多光谱遥感传感器，具有36个观测通道，覆盖了当前主要遥感卫星的主要观测数据。利用MODIS的陆地产品提取秸秆焚烧火点的处理流程如下：
图2：技术流程图
流程涉及到影像的光谱处理和几何处理。ENVI/IDL遥感影像处理软件能对MODIS数据进行快速精确的光谱和几何处理，给更快速、更准确的监测秸秆焚烧现状提供了技术保障。其中大气校正模块FLAASH是采用精度最高大气校正算法，是基于光谱信息提取相关专题产品的必要前提；其面向对象信息提取工具ENVI EX能有效的对地表进行分类，提取土地覆盖信息；其强大的波段运算工具是提取各种专题数据的有效手段。
由于MODIS的 MOD14热异常数据可供直接获取使用，能够探测比气象卫星更小更多的火点（面积50平方米），是监测秸秆焚烧理想的数据源。下面从数据源获取、火点提取、农用地提取、结果输出等几个方面来介绍ENVI在秸秆焚烧监测中的应用方案。
1、数据获取
MODIS数据下载地址：http://ladsweb.nascom.nasa.gov/data/search.html。
(一) 首先打开网页，选中需要的数据源类型。其中，日期类型为：月/日/年 时:分:秒，其中，网页中显示的时间为UTC时间。
(二) 在“spatial selection” 选项中选择“latitude/longtitude”，按监测区域的经纬度选择影像范围。点击“search”查到需要的数据，勾选需要的数据，点击“order files now”， 输入你接收信息的邮箱，点“order”开始订购该数据。(如果要搜索多天数据，可以选“add files to shopping cart”继续搜索其他日期的数据。
(三) 所订购数据的存放位置信息： 点击“Data->Track Orders ”可以查看所有已订购的数据的状态。如果“state”显示“avalable”即可开始下载。
本实验数据覆盖范围：
图3：数据覆盖范围示意图
2、火点提取
图4：Modis数据
说明：在这个步骤中，一般情况下，如果不是采用MOD14数据，而是由MOD02数据以及MOD03数据进行热异常信息提取，则需要用到ENVI的FLAASH大气校正模块进行大气校正。在这里，我们简化地采用MOD14数据。
第一步：MOD14数据的地理定位
MOD14热异常数据地理定位必须要用到MOD03数据。
1、 “MOD14.A2012163.0305.005.2012163094359.hdf”显示如下：
图5：MOD14热异常数据集
2、 选择MOD03的经纬度数据打开：
图6：MOD03数据集信息
3、利用经纬度构建GLT地理查找表
图7：GLT文件构建
“X Geometry Band”即为“Longitude”， “Y Geometry Band”即为“Latitude”：选择输出坐标系,设置GLT参数，创建完成后，在Available Band List出现。
图8：GLT文件构建
4、利用GLT地理查找表对MOD14数据进行地理定位。
利用上一步生成的GLT地理查找表对MOD14进行地理定位，坐标系为WGS-84，结果如下：
图9：地理定位后的MOD14热异常数据
第二步：火点提取
对MOD14地理定位后的波段进行波段信息统计，统计信息如下：
图10：MOD14统计结果
0-9分别代表的含义如下表所示：
图11：MOD14数据像元值分类
可见，像素值为7，8，9的地方为火点。利用波段运算窗口进行火点提取：
BandMath表达式：(b1 ge 7)*0+(b1 lt 7)*1
注：就是讲7、8、9像元值变成1，其余为0，对MOD14数据进行二值化处理。
图12：提取火点信息的BandMath表达式
在表达式窗口中输入二值化表达式对图像进行二值化，火点区赋值为0，非火点区赋值为1. 得到的结果如下图所示：
图13：MOD14二值化结果
利用感兴趣区域将火点影像转化为矢量，结果（红色三角形点）如下：
图14：疑似火点矢量化结果
3、农用地提取
说明：在此步骤中，一般情况下，需要使用实时 MOD02数据进行农用地提取，此时需要用到ENVI EX面向对象信息提取工具对土地覆盖进行分类提取。在这里我们简华地采用现势性不强（2009年）的MCD12Q1数据作为替代。
火点提取出来后，接下来需要判定火点是否在农用地之上，在此，我们认为在农用地之上的火点为秸秆焚烧造成的火点。农用地的提取方法有多种，在此我们重点介绍由MODIS三级数据土地覆盖类型产品提取农用地的方法，其余方法则不一一列举。
MODIS三级数据土地覆盖类型产品（Land Cover data）是根据一年的Terra和Aqua观测所得的数据经过处理，描述土地覆盖的类型。该土地覆盖数据集中包含了17个主要土地覆盖类型，根据国际地圈生物圈计划（IGBP），其中包括11个自然植被类型，3个土地开发和镶嵌的地类和3个非草木土地类型定义类。
MODIS Terra + Aqua三级土地覆盖类型年度全球500米产品MCD12Q1采用五种不同的土地覆盖分类方案，信息提取主要技术是监督决策树分类。下面是该数据中包含的五个数据集，五个分类方案如下：
土地覆盖分类1：IGBP的全球植被分类方案
土地覆盖分类2：美国马里兰大学（UMD格式）方案
土地覆盖分类3：基于MODIS叶面积指数/光合有效辐射方案
土地覆盖分类4：基于MODIS衍生净初级生产力（NPP）方案
土地覆盖分类5：植物功能型（肺功能）方案
MCD12Q1是HDF格式的，包含16个数据集
图15：MCD12Q1数据集说明
第一步：数据读取
我们利用MCD12的第一类分类数据进行农用地信息提取。
图16：MCD12Q1数据集信息
（1）Land_cover_type数据说明
打开Land_cover_type1到Land_cover_type5，是五类分类方案的分类图，图像上的DN值各个值都代表某一地物类型，具体意义如下表：
图17：MCD12Q1数据分类说明
我们感兴趣的是第一种分类图像中DN值为12的像元。
第二步：数据地理定位
MODIS三级数据产品土地覆盖数据是已经经过投影了的，其投影类型为正弦曲线投影，但是很多遥感软件在打开HDF数据的时候，会有投影丢失现象，ENVI也是，可以通过hdf属性查看工具，查看ENVI打开hdf时未识别的投影信息，以及投影参数，利用该参数在ENVI中自定义正弦曲线投影，为MCD12Q1数据添加投影。
（1）用hdf属性查看工具查看该hdf数据的地理定位信息，可得图像左上角、右下角坐标值，投影类型，投影参数，地球半径，起始位置。
图18：MCD12Q1数据集的地理定位信息
（2）根据投影信息，在ENVI中自定义投影：
图19：自定义正弦曲线投影坐标
并修改其头文件信息：
图20：编辑MapInfo坐标信息MCD12Q1数据
已经具有投影信息的数据：
图21：定义Mapinfo坐标信息后的MCD12Q1数据
第三步：投影变换与镶嵌
（1） 针对每一张MODIS土地覆盖的影像，经过地理定位后，将其坐标转到WGS-84坐标系下，用投影转换工具对相应的影像进行投影转换，变换到WGS84坐标系下：
图22：MCD12Q1数据投影转换
（2）影像镶嵌
所有的土地覆盖影像进行投影转换之后，将此数十幅影像进行拼接。用基于地理投影的镶嵌工具进行影像镶嵌，如下图所示：
图23：MCD12Q1数据镶嵌面板
将所有的背景值设为0，避免互相压盖。镶嵌结果如下：
图24：MCD12Q1数据镶嵌结果
同样，利用火点提取一节中的方法，依照MCDQ1的分类说明，把DN=12的像素（农耕用地）提取出来。结果如下图所示（绿色区域）
图25：从MCD12Q1数据集中提取农耕用地
4、结果输出
在Image窗口中将结果连接至ArcMap，第二个环节的火点与第三环节的土地覆盖叠加显示如下：
图26：火点提取结果
该结果可以在ArcGIS中制图制图输出或发布服务。
5、IDL定制开发
ENVI是基于IDL语言开发的，也提供了丰富的开发接口供用户调用，上述流程在ENVI中都有相应的API，利用IDL语言将各个流程串接起来，可以形成专门的工具嵌入ENVI主界面，或作为一个独立的IDL程序提供给软件使用者，亦可嵌入其他语言做系统集成开发，可以降低对遥感专业知识的要求，直接从原始数据到最终结果。当然，开发需要专业的开发人员参与。
下面是基于ENVI/IDL开发的基于HJ的火点检测的程序（2009中南大学，彭光雄），火点检测采用如下算法（结果与MOD14相当）：
此时，只需要点击菜单，选择影像：
输入少数参数：
就可以得到火点分布情况：
该程序可以作为ENVI的一个插件，也可以作为一个单独的应用程序，或是与C# /Java做混合开发嵌入业务系统中。
文章来源：http://blog.sina.com.cn/s/blog_764b1e9d01010sl9.html
