基于ENVI和哨兵2数据提取云南玉溪和安宁山火受灾面积
引言
2023年4月11日15时27分，云南省玉溪市江川区九溪镇发生森林火情。当地森林草原防灭火指挥部立即启动应急预案，组织力量扑救。省、市、区各级共投入森林消防、消防救援、公安民警、武警官兵、专业扑火队、民兵、干部群众共计4000余人开展扑救工作，妥善转移安置部分可能受火灾影响村庄的群众。经过全体扑救人员100多个小时艰苦鏖战，15日21时55分，玉溪市江川区“4·11”森林火灾明火已扑灭，现已转入火场清理看守阶段。此次森林火灾无人员伤亡。
图1：4月13日云南省消防救援总队的消防指战员在玉溪市江川区火场扑救山火。新华社发
2023年4月13日16时21分，云南省安宁市温泉街道官庄村发生森林火情，安宁市第一时间启动应急预案，成立前线指挥部，云南省、昆明市有关部门赶赴现场指挥扑救工作。目前，云南省、昆明市、安宁市共投入森林消防、消防救援、民兵、干部群众等1400余人，出动直升机40架次，工程机械、水车等56台（辆）开展扑救。云南省安宁市森林草原防灭火指挥部消息，截至4月18日8点26分，昆明市安宁市“4·13”温泉街道官庄村森林火灾明火已扑灭，全部力量转入火场清理看守，无人员伤亡情况。
图2：森林消防指战员在安宁火场进行灭火作战。张瑞坤摄
通过卫星遥感手段，可以监测火灾、预测火灾的发展趋势。获取2023年4月16日的Landsat8影像，通过短波红外的假彩色显示方法，可查看到安宁山火当时的火点区域，如下图所示：
图3：采用R:SWIR(2200nm)、G:NIR(865nm)、B:BLUE(483nm)颜色加载方案，橘红色为火点区域
1 数据源
对于灾后林火受灾面积提取，可以通过前后两时期影像，使用ENVI直接变化检测工作流工具进行提取。
注：可访问 envi.geoscene.cn/envi_license 获取最新ENVI5.6.3软件试用。中文语言包可在app store中下载安装，App Store 官网：envi.geoscene.cn/appstore。
数据源采用2023年4月11日和4月21日山火前后，两期哨兵二号L2级的影像。哨兵二号L2级的影像已经进行几何校正和辐射校正处理，为地表真实反射率数据。
- S2B_MSIL2A_20230411T033539_N0509_R061_T47RRH_20230411T065637
- S2B_MSIL2A_20230421T033539_N0509_R061_T47RRH_20230421T063920
以上数据可从百度网盘下载，下载链接：
链接：https://pan.baidu.com/s/1s0wv_6goLF-rzGXTIV0yPw?pwd=envi
提取码：envi
2 数据预处理
2.1 打开数据
在ENVI菜单栏，选择“文件>打开…”定位到解压后的哨兵二号影像文件夹（可解压到一个短目录），分别选择文件夹内的“MTD_MSIL2A.xml”元文件打开两期数据。在左侧Layer Manager图层管理分别选中哨兵二号影像，右键选择“修改RGB波段”，依次选择B8近红波段、B4红波段、B3绿波段，使用标准假彩色显示方法，植被会显示红色。观察两时期影像，可以看到火灾受灾区。
图4：火灾受灾区域
2.2 研究区裁剪
研究区裁剪首先需要绘制感兴趣区域。在ENVI上方工具栏，选择ROI感兴趣工具，图形选择矩形，分别绘制火灾区域矩形范围。
图5：绘制研究区范围
绘制好ROI文件之后就可以进行矩形范围裁剪。在ENVI菜单栏，选择“文件>另存为>另存为ENVI”，选择10m-S2MSI2A 2023-04-11…哨兵二号10米4波段数据集，之后选择空间裁剪，选择ROI文件，选择玉溪山火矩形范围文件进行矩形范围裁剪。同样的方法处理另一时期的影像和安宁山火区域对应的两期影像。
图6：裁剪影像
裁剪后得到安宁市和玉溪市山火矩形范围影像，如下图所示：
图7：安宁火灾区域（上）和玉溪火灾区域（下）
2.3 短波红外波段组合（选做）
短波红外组合主要是为了在进行图像直接变化检测的时候可以构建燃烧指数(Normalized Burn Ratio，NBR)进行火灾范围提取，构建此指数需要用到短波红外波段。也可以跳过此步，直接使用NDVI指数进行火灾范围提取。在ENVI右侧Toolbox工具箱，选择栅格数据管理>波段组合工具。在参数设置对话框，首先选择上一步裁剪得到的玉溪山火地区10米4波段数据，之后选择哨兵2号20米6波段数据，其他参数按照默认，设置输出路径和文件名，点击确定按钮。得到包含短波红外波段的波段组合结果。
图8：短波红外波段组合
3 图像直接比较法变化检测
直接变化检测工作流（Image Change Workflow）工具可以检测两个时相影像中增加和减少两种变化信息，适合获取地表绝对变化信息。利用直接变化检测工作流工具提取火灾范围。
第一步、启动直接变化检测工作流
在工具箱列表中，双击/变化检测/直接变化检测工作流，打开文件选择对话框，分别为前一时相图像选择S2B_MSIL2A_20230411_玉溪_短波红外波段组合.dat和后一时相图像选择S2B_MSIL2A_20230421_玉溪_短波红外波段组合.dat。
图9：直接变化检测工作流工具
单击下一步按钮打开图像配准面板。直接单击下一步，跳过图像配准步骤。
第二步、变化信息检测
在检测方法选择面板中，提供两种方法：图像差值法和图像变换法。选择图像差值方法单击下一步。
图像差值法又提供了三种方法：
- 波段差值
在选择输入波段选择相应的波段。切换到高级，提供辐射归一化（Radiometric Normalization）选项，可以将两个图像近似在一个天气条件下成像（以Time1图像为基准）。
- 特征指数差
这个方法要求数据是多光谱或者高光谱，自动根据图像信息（波段数和中心波长信息）在选择特征指数列表中可以选择的特征指数。提供四种特征指数：
- 归一化植被指数Vegetation Index (NDVI)：植被区域NDVI值大。
- 归一化水指数Water Index (NDWI)：水体区域NDWI值大。
- 归一化建筑物指数Built-up Index (NDBI)：建筑物区域NDBI值大。
- 燃烧指数Burn Index：燃烧区域BI值大。
图10：变化检测方法选择面板
切换到高级，自动为波段1和波段2选择相应的波段（前提是有中心波长信息），否则手动选择。
- 选择特征指数，在选择特征指数列表中选择植被指数 (NDVI)，如果有短波红外波段也可选择燃烧指数。
- 勾选预览选项，可以预览变化信息检测的结果，红色区域表示为NDVI减少，蓝色区域表示NDVI增加。
- 波谱角
光谱角度差异通过光谱角度进行差异分析。该选项通过计算两期影像像元光谱曲线之间的角度来确定前一时期影像光谱和后一时期影像光谱之间的相似性。对于高光谱数据，推荐使用光谱角度差法。如果单个输入波段没有足够的信息，或者对于波段差值法来说有太多的波段可供选择，或者对于特征指数差法来说没有感兴趣的特征，那么这是一个不错的选择。
这里选择特征指数差法，单击下一步打开阈值与输出面板，提供两种方法：
- 应用阈值：设置阈值获得增加减少区域结果。
- 仅输出差值图像：直接输出差值图像。
选择应用阈值，单击下一步打开修改阈值面板。
第三步、变化信息提取
直接变化检测工作流工具可以从变化检测结果中可以提取三种变化信息：
- Increase and Decrease：增加 (蓝色) 和减少(红色)变化信息
- Increase Only：增加 (蓝色)变化信息
- Decrease Only：减少(红色)变化信息
提供两种阈值设置方法：
- 自动设置阈值（Auto-Thresholding）
在选择感兴趣的变化列表中选择仅减少，获取NDVI减少区域。
提供四种算法自动获取分割阈值：
- Otsu's: 基于直方图形状的方法。使用直方图积累区间来划分阈值。
- Tsai's: 基于力矩的方法。
- Kapur's: 基于信息熵的方法。
- Kittler's: 基于直方图形状的方法。把直方图近似高斯双峰从而找到拐点。
切换到手动设置阈值查看获取的分类阈值，可以手动更改，勾选预览查看分类效果。
在选择自动阈值方法列表中选择Otsu's。
图11：修改阈值面板
- 手动设置阈值（Manual）
在手动选项卡，可以根据自动设置的阈值上下微调，获得更为准确的结果。此处设置为09。
勾选预览查看效果。
单击下一步打开小斑清除优化结果面板。这个面板的作用是移除椒盐噪声和去除小面积斑块。
- 勾选平滑，平滑核（Smooth Kernel Size）：5，值越大，平滑尺度越大。
- 勾选聚合，最小聚类值（Aggregate Minimum Size）：100。
勾选预览查看效果，单击下一步打开输出保存结果面板。
第四步、输出变化信息
可以输出四种结果或格式：
- 以图像格式输出变化结果
- 以矢量格式输出变化结果
- 变化统计文本文件
- 输出差值图像
勾选输出变化结果矢量文件，选择输出为Shapefile格式。
切换到其他输出，勾选输出变化结果统计值，输出统计文件。
单击完成输出结果。
图12：变化检测工作流提取的火灾区域
第五步、火灾区域矢量导出
如果变化检测工作流提取的有多余的变化地物，可以选择想要的变化图斑重新导出为一个新的矢量文件。按住键盘的Ctrl键，鼠标左键选中要输出的矢量图斑，在左侧的图层管理选中变化检测得到的矢量结果，右键>查看/编辑属性。
图13：编辑属性
在属性浏览器，选择文件>保存选中记录为新的Shapefile。
图14：保存选中矢量数据
得到最终的火灾受灾范围如下图所示：
图15：玉溪林火受灾范围
通过矢量属性里面的面积字段统计得到玉溪山火火灾受灾面积为：6124.97 公顷。制图结果如下图所示：
图16：玉溪林火受灾范围制图
使用NDVI差值、燃烧指数差值和波谱角度差三种方法提取的玉溪火灾受灾区如下图所示：
图17：三种方法提取火灾受灾面积
同样的方法提取安宁山火受灾范围，最终两地山火受灾范围结果如下图所示：
图18：云南安宁山火和玉溪山火受灾范围提取结果
4 总结
ENVI图像直接变化检测工作流工具可以方便快速地提取山火受灾范围，为灾后评估和生态恢复提供数据支撑。此外，在山火救灾过程中实时的火点数据则可以提供灭火指引，为消防人员提供火灾的位置和规模信息，为火灾救灾提供帮助。
