长篇幅详解辐射定标、大气校正、监督分类、掩膜统计、植被覆盖度操作

一、序言

Landsat卫星是遥感初学者甚至遥感从业者、相关研究学者接触最多的一个遥感数据源。今天我们通过一个示例来学习ENVI一些相关的操作。

1.1 实验区简介

大家在浏览谷歌地球或者谷歌地图的时候，有没有发现陕西北部的卫星图会有上图这样的一个明显的轮廓特征。这个轮廓鲜明的地区就是陕西的吴起县，如图下所示：

1998年，对吴起来说是一个至关重要的年份：这一年，该县率先在全国开启退耕还林探索，在全县范围内推行“封山禁牧、舍饲养畜”政策，一场旨在改变生存环境的生态环境建设在吴起县轰轰烈烈展开。20年，吴起累计退耕还林243.79万亩，因封得最早、退得最快、面积最大，成为全国退耕还林第一县。这个卫星图上清晰可见的轮廓正是人们改造自然的见证。

1.2 数据查询选取

首先打开ArcMAP，加载全国的县级区划矢量，使用属性查找功能选择吴起县：

这样吴起县便被选中，然后右键可以导出选中的吴起县的区划为单独一个SHP文件：


然后加载Landsat的PATH和ROW的矢量文件，标注出PATH和ROW：

我们可以知道覆盖吴起县需要128034和128035两景Landsat数据：

接下来我们可以将当前视图的图框坐标系设置为UTM-48N：

打开编辑，画一个比吴起县大的矩形，作为我们的研究区范围（因为我们需要吴起县和县域外一部分数据做比对，所以就没有直接用吴起县的区划矢量），然后导出这个矩形为一个新的SHP文件：

接下来可以进入USGS查询数据，USGS可以导入研究区的SHP文件，将SHP文件压缩为ZIP文件即可被导入：


接下来选择卫星查询即可，我们本次选用的数据是1997年5月19日的Landsat5和2018年6月14日Landsat8。USGS查询数据的时候，可以先通过贴图的方式来看看选取的数据是否理想，比如云量这些：

二、预处理

下面我们以2018年的两景数据为例来介绍一般的预处理过程，涉及数据裁切、辐射定标、掩膜的使用、大气校正等操作。首先打开ENVI软件打开数据，常规的方法是使用OpenAs里面找到Landsat卫星：

偷懒的方法是将解压后文件夹里的MTL.txt文件直接拖拽进ENVI视窗即可：

切忌自己使用Laystack工具私下将这些波段组合在一起，不仅费事费力，而且自己组合后的文件不带定标参数，做不了后续操作。

2.1辐射定标

常规的方法一般是对整景数据进行定标，但是这样的话，假如研究区比较小，反而会增加数据量和运行时间，所以我们可以先裁剪，再进行辐射定标。
裁剪的时候需要将MASK选项设为YES，这样就是根据实际形状裁剪，默认的NO则是裁剪外接矩形

裁剪完的数据，定标参数都在，我们可以放心的进行辐射定标了：

下面对裁剪的后的数据进行辐射定标的操作，记得选择Aplly FLAASH Settings，输出路径尽量不要中文：

接下来我们对上面另一景数据也进行辐射定标操作，其实在辐射定标的工具里面有一个裁剪的选项，我们可以连带将裁剪和辐射定标一起做了：

2.2掩膜统计

辐射定标后一般是进行FLAASH大气校正，因为大气校正的选项里有一个平均高程的选项，所以我们穿插进来一个掩膜统计的步骤。在ENVI的安装目录里自带一个全球的DEM数据，分辨率为900米，不想下载DEM的话，可以使用这个进行掩膜统计平均高程。
首先打开这个数据后，点击工具箱里的统计工具，在MASK option里选择Build MASK，定义MASK的时候可以选择导入我们自己画的ROI如下图：

也可以选择我们的研究区范围的SHP文件：



导入ROI或者SHP后就可以进行掩膜统计了：

下图是进行掩膜统计的结果，平均高程大概是1500米，其实与百度的结果差不多：

2.3大气校正

对于不同拍摄时间，不同条带的数据需要分别辐射定标过后单独做大气校正，对于同一条带同一天拍摄的上下两景数据，因为一般两景的拍摄时间间隔很短，大气条件几乎一样，可以先拼接再进行大气校正。

2.3.1 单景大气校正
我们首先来介绍单景的大气校正，FLAASH大气校正的输入输出路径尽量简洁，短路径不带中文，平均高程的话单位是KM，刚才我们统计的高程是1500米，所以这里填1.5KM，飞行时间记事本或者写字板打开MTL文件里面有景中心时间，另外KT那里别忘了选：


分块处理也最好关掉，不容易报错：

然后对另一景采用同样的操作。

大气校正后查看光谱曲线，我们可以插入一个新视图，视图1打开的是大气校正前的数据，视图2打开大气校正后的数据：

先放大数据找到一块植被区，点击光标工具，复制GEO坐标：

复制后的GEO坐标粘贴进去搜索窗口，点击光谱曲线按钮，既显示出此坐标点的光谱曲线：

对另一个视图进行相同的操作，这样大气校正前后的光谱曲线就分别显示出来了：

2.3.2 拼接后大气校正

这次涉及的数据是同一条带同一天拍摄的上下两景数据，我们可以先将辐射定标后的数据进行拼接再大气校正。由于本区域上下两景数据跨带了，一个是48一个是49，在拼接前需要先进行重投影，本次我们选择48带，因此需要将128034这个数据重投影为UTM 48N：

打开重投影工具，输出坐标系设置如下：

重投影后辐射定标的数据类型由原来的BIL变为BSQ，因此我们还需要对其进行数据类型转换，因为FLAASH大气校正的输入数据类型要求为BIL：

对于两个数据叠加出现黑色背景的情况，比如下图，需要先进行编辑头文件，添加忽略背景值：

对于编辑头文件出现报错的情况，可以进classic模式，这样就不会报错了：

忽略背景值后两个数据正常叠加了：

启动镶嵌工具进行拼接，可以生成一个拼接线：

编辑拼接线，尽量沿着线性地物画拼接线：

同一天拍摄的数据不选直方图匹配：

镶嵌后的结果：

接下来进行大气校正，参数的填写与上一篇类似，只是飞行时间那里，哪一景占的比例比较大就填写那一景的飞行时间。

三、监督分类

数据镶嵌后就可以进行监督分类了，对于目视选取样本，我们需要选择对于目视判读最合适的组合方式，下图是本区域Landsat8数据不同的波段组合显示效果：

本次我们选用654组合选择样本，当然也可以随时切换其他组合一起判读：

结合实地的自然条件，我们选择有林地、坡耕地、水浇地、草地、水域、居民地、裸地七个类别进行监督分类。
本次数据拍摄时间是6月14日，林地长势茂盛在654组合下显示为鲜绿色：

由于陕北7、8月才进入雨季，坡耕地上基本无有效灌溉，因此6月份坡耕地的作物长势还还很弱，显示为红色：

地势平坦具备灌溉条件的耕地我们命名为水浇地，呈现鲜明的耕地纹理特征，颜色鲜绿：

可能大家做监督分类最头痛的是草地样本的选取，大城市好一些还有一些公园球场之类的，或者西北、西南、青藏这些典型的草原区，但是在其他草地分布并不突出的地区，草地确实比较难选，草地会零星分布在一些河滩，河边这些地方，草地没有耕地那么典型的纹理：


水域样本：

居民地样本：

裸地样本：

本次每一个类别的样本均选取了15个以上，保证分布均匀，然后计算可分离度，可分离度1.8以上即为良好：

采用SVM分类方法：

四、植被覆盖度反演

这一章我们讲植被覆盖度反演操作，需要用到NDVI和监督分类的结果。

4.1 计算NDVI

计算NDVI有两种方法，一种是使用BandMath，一种是使用NDVI工具，不管是哪种方法都要选对波段：



这两种方法计算出来的结果是一样的：

计算好的NDVI需要统计一下有没有异常值，假如存在异常值可以波段运算中输入公式：b1>-1<1,去除异常值，不规则研究区统计可参考前文统计DEM的方法，只对研究区的范围进行掩膜统计，排除背景的影响。

4.2 构建掩膜

这一步我们需要将用先前监督分类的结果，将每一个地类制作为掩膜文件。首先打开ApplyMask，选择NDVI为基准图像，选择BuildMask选项，首先构建有林地掩膜：

导入掩膜数据的时候选择监督分类的结果：

有林地的值为1，所以我们这里值范围都为1：


接下来依次构建，坡耕地、水浇地、草地、居民地、裸地的掩膜，操作与上面类似，只是值域范围需要特别注意，不要搞混：

4.3 NDVI最大值最小值统计

接下来使用统计工具统计每个地类的NDVI区间，首先统计有林地的NDVI范围：



确定最大值最小值可以选择直方图的拐点处，也可以选择2%-98%水平的置信区间，也可以选择统计值的突变处，比如统计值突然由三位数变成四位数或者由二位数突变为三位数的值处：


下图是每个地类本次的统计结果

4.4 植被覆盖度计算

首先计算NDVIsoil，输入公式：b10.5406+b20.1914+b30.4151+b40.4638+b60.0331+b70.077


然后计算NDVIveg，一样的方法，输入公式：b10.8912+b20.6951+b30.6992+b40.7124+b60.5713+b70.4096

最后计算植被覆盖度

对植被覆盖度结果进行去除异常值处理：b1>0<1，这样可以去除异常值和水体中的Inf：

五、数据展示

上面讲了我们对研究区2018年的数据进行的操作流程，同样的方法对1997年数据进行相同的操作，得到两期的数据成果。
我们对两期的监督分类数据分别单独将林地提取出来，波段运算输入b1 eq 1，将林地提取出：

分别使用New color slice选项，对提取出的林地进行颜色渲染，叠加行政区，可以明显看出20年来，吴起县的林地扩张情况：
也可以使用掩膜统计进行面积统计，本次就不再说明了。

同样也可以导入GIS制图，我们以植被覆盖度为例导入GIS操作：

首先，ENVI里File下Save as TIF文件，然后使用Arcmap加载，
符号系统里点击分类，设为统一的间断：

对于初次打开TIF文件，提示计算直方图一定要选YES，这样才能正常显示图像：

选择色带，加入标注：

切换布局视图，插入数据框，插入图例、指北针等：

导出地图即可：

好了今天介绍就到这里，有兴趣的可以关注下，博文主页：
https://blog.csdn.net/qq_46071146
更多好文章：