ENVI基于Landsat影像提取水体流程

一、导入数据

启动ENVI,file—Open As,根据影像选择合适的传感器，此处以Landsat8为例，因而选择光学传感器下Landsat带有源数据的tiff图；File—Open As—Optical Sensors—Landsat—Geo TIFF with Metadata,导入MTL.txt文件(影像元数据，头文件)

※※※ 数据如果遇到打不开的情况，查看MTL.txt文件中的开头和结尾：

※ 对于L1GS类型数据：

…………………………………

若MTL文件中开头为：

GROUP = LANDSAT_METADATA_FILE

GROUP = PRODUCT_CONTENTS

……………………………………

MTL文件中结尾为：

……………………………………

END_GROUP = PRODUCT_PARAMETERS

END_GROUP = LANDSAT_METADATA_FILE

END

………………………………………

则修改GROUP = LANDSAT_METADATA_FILE为GROUP = L1_METADATA_FILE

END_GROUP = LANDSAT_METADATA_FILE为 END_GROUP = L1_METADATA_FILE

修改后保存txt，即可打开。

※ 对于L1GS类型数据：

修改GROUP = LANDSAT_METADATA_FILE为GROUP = L1_METADATA_FILE

END_GROUP = LANDSAT_METADATA_FILE为 END_GROUP = L1_METADATA_FILE

另外删掉GROUP = LEVEL1_PROCESSING_RECORD

到END_GROUP = LEVEL1_PROCESSING_RECORD之间的数据（若删掉后仍然打不开，则保留GROUP = LEVEL1_PROCESSING_RECORD和END_GROUP ==LEVEL1_PROCESSING_RECORD这两句话，只删掉这两句话中间的内容）

若上述方法仍旧打开失败，则考虑升级ENVI版本，例如升级至ENVI5.6；

若以上方法均无法导入MTL.txt，则需要进行波段合成。

二、图像预处理

1.几何校正

（Landsat7和8不需要进行几何校正）

几何校正是利用地面控制点和几何校正数学模型来矫正非系统因素产生的误差，校正过程中会将坐标系统赋予图像数据。ENVI的几何校正功能中包括控制点选择方式和几何校正模型。控制点选择方式包括：

• 从栅格图像上选择：可以从经过校正的影像、地形图等栅格数据中选择控制点，对应的控制点选择模式为Image to Image；

• 从矢量数据中选择：可以从经过校正的影像、地形图等矢量数据中选择控制点，对应的控制点选择模式为Image to Map；

• 从文本文件中导入控制点：可以从事先经过GPS测量、摄影测量或者其他途径获得的控制点的坐标数据生成的[Map (x,y), Image (x,y)]格式的文本文件直接导入作为控制点，对应的控制点选择模式为Image to Image 和Image to Map；

• 键盘输入：通过键盘敲入坐标数据并在影像上找到对应点。

        1）从栅格图像上选择

        2）从矢量数据中选择

        3）从文本文件中导入

        4）键盘输入

几何校正模型包括：

• 仿射变换（RST）：包括平移（Translation）、旋转（Rotation）、缩放（Scaling）、剪切（Shearing）等；

• 多项式模型：应用此模型时，需要确定多项式的次方数，通常选择2次或3次。

• 局部三角网（Delaunay Triangulation）：

操作：控制面板File—open，打开需要校正的图像；工具箱中，选择Geometric Correction—Registration—Image Registration Workflow，启动几何校正模块。

在Base Image File 对话框中，选择基准图像(已校正好的图像)；在Warp Image File 对话框中，选择待校正图像，点击Next ；

1）生成控制点：

①Main 选项卡参数设置

Matching Method（匹配算法）：​Cross Correlation 方法一般用于相同形态的图像，如两幅光学图像；​Mutual Information 方法一般用于不同形态的图像，如光学-雷达图像；

Minimum Matching Score（最小匹配点匹配阈值）：ENVI 的自动找点功能找到的同名点​低于这个阈值，则会自动删除而不参与校正，阈值范围是0～1；

Geometric Model（几何模型）：过滤匹配点的几何模型；FittingGlobal Transform 方法适合于绝大部分图像，Frame Central Projection 方法适合框幅式中心投影的航空图像数据； 

Transform（变换模型）：包括一次多项式（First-Order Polynomial）和仿射变化（RST）；

Maximum Allowable Error Per Tie Point（每个匹配点最大容许误差）：值越大，精度越差。

②Seed Tie Points 选项卡参数设置

点击Start Editing 按钮手动添加控制点，此时鼠标变成+号；

添加和删除控制点：

点击“Add”按钮，在一副影像上选择好特征点位后单击左键—右键选择“Accept as Individual Points—跳转至另外一幅影像，在这幅影像上添加相应点，至此控制点添加完成，重复上述工作添加多个控制点；

点击“Delete”按钮，选中要删除的控制点，右击删除；

查看控制点：停止编辑—点击Show Table —查看控制点

③Advanced 选项卡参数设置

Matching Band in Base Image：基准图像配准波段；

Matching Band in Warp Image：待配准图像配准波段；

Requested Number of Tie Points：拟匹配的连接点个数；

Search Window Size：搜索窗口的大小，需要大于匹配窗口；

Matching Window Size：匹配窗口的大小；

Insert Operator：匹配算子，forstner 方法精度较高，速度最慢。

检查控制点和预览结果：

Tie Points模块中打开Show Table，其中IMAGE1X、IMAGE1Y 为基准图像上连接点的X 像素坐标（列数）、Y 像素坐标（行数）；IMAGE2X、IMAGE2Y 为待校正图像上连接点的X 像素坐标（列数）、Y像素坐标（行数）；ERROR 为误差值，选中ERROR 列右键选择Sort by Selected Column Reverse按照误差从大到小排序，找出误差大于1 的点删除。

2.地形校正

（Landsat7和8不需要进行地形校正）

ENVI App Store中下载ENVI地形校正工具：

下载地址：

https://envi.geoscene.cn/appstore/toprhcrion/

安装方法：

将压缩包内的所有文件夹（如 extensions、custom_code 等）拷贝到

• ENVI 5.3：C:\Program Files\Exelis\ENVI53\
• ENVI 5.4~5.7：C:\Program Files\Harris\ENVI5x\
• ENVI 6.0~：C:\Program Files\NV5\ENVI6x\

下的ENVI 安装路径中，覆盖同名文件夹，重启ENVI。

加载需要校正的图像：

File—Open

地形校正：

Toolbox—Extensions—Topographic Correction，在地形校正控制面板中输入以下参数：

Input Spectral Raster：输入待地形校正影像；

Input DEM Raster：输入DEM影像；该工具会自动按照待校正影像对DEM影像进行重投影和重采样，使其与待校正影像空间分辨率和投影一致；

Sun Azimuth：太阳方位角，元数据（******_MTL.txt文件）可查；

Sun Elevation：太阳高度角，元数据（******_MTL.txt文件）可查。

Kernel Size：卷积核大小。地形校正需由DEM数据计算坡度、坡向，卷积核的值越大，计算出的地形因子越平滑；

Grid Size：参与回归运算样本点采样间隔，默认为10，即每10个像元取一个样本点；若需全部像元参与回归，则设置为1；

Topo Model：地形校正模型；

Resample Method：DEM数据重采样方法；

Output Raster：输出校正结果；

Display Result：在ENVI中显示校正结果。

3.辐射定标

目的：将记录的原始像元亮度DN值（记录地物灰度值，即传感器得到的数字测量值）转换为地表反射率或与表面温度等物理量有关的相对值；将图像的亮度灰度值转换为绝对辐射亮度。

操作：找到envi界面右侧的Toolbox–Radiometric Correction(辐射校正)–Radiometric Calibration(辐射定标)

选择多光谱图像进行定标，点击OK

4.大气校正

Toolbox–Radiometric Correction(辐射校正)—Atmospheric Correction Module(大气校正模块)—FLAASH(大气校正（FLAASH Atmospheric Correction），打开大气校正参数面板

• Input Radiance Image：选择上一步辐射定标后的图像，在弹出的对话框选里面选择use single scale factor for all bands；

• Output Reflectance File：输出反射范围文件，设置输出路径；

PS:输出结果默认是经过缩放的反射率数据（放大了10000倍），其像素值的范围从0到10000（表示0到100%的反射率）。一些像素可能位于该范围之外，这些异常值通常位于饱和度过高的高反射表面，或者暗像元。负值通常出现在反射率低的深水或阴影中。要将像素值缩放为0到1之间，可使用Band Math工具进行波段运算，公式为b1/10000.0，将b1指定为整个文件。

• Output Directory for FLAASH Files：大气校正文件的输出目录；

• Rootname for FLAASH Files：大气校正文件命名；

• Sensor Type：传感器类型，选择影像对应的类型；

• Flight Date：时间，右键图层——view metadata——time可查看时间；

• Ground Elevation：平均海拔，查看平均海拔：可以通过百度搜索或参考如下步骤：

         File-open,找到ENVI安装路径下ENVI自带的全球900m DEM数据文件，打开

 Tooxbox中搜索Compute Statistics，双击打开

        进行如下操作：

依次点击Subset by File Input File、Select Statistic Subset和Compute Statistics Input File中的OK

计算统计面板参数设置如下：

统计结果如下：

• Atmospheric Model：大气模型，根据下表选择合适的大气模型

Latitude	1月	3月	5月	7月	9月	11月
80°N	SAW	SAW	SAW	SAW	SAW	SAW
70°N	SAW	SAW	MLW	MLW	MLW	SAW
60°N	MLW	MLW	MLW	SAS	SAS	MLW
50°N	MLW	MLW	SAS	SAS	SAS	SAS
40°N	SAS	SAS	SAS	MLS	MLS	SAS
30°N	MLS	MLS	MLS	T	T	MLS
20°N	T	T	T	T	T	T
10°N	T	T	T	T	T	T
0°	T	T	T	T	T	T
10°S	T	T	T	T	T	T
20°S	T	T	T	MLS	MLS	T
30°S	MLS	MLS	MLS	MLS	MLS	MLS
40°S	SAS	SAS	SAS	SAS	SAS	SAS
50°S	SAS	SAS	SAS	MLW	MLW	SAS
60°S	MLW	MLW	MLW	MLW	MLW	MLW
70°S	MLW	MLW	MLW	MLW	MLW	MLW
80°S	MLW	MLW	MLW	MLW	SAW	MLW

• aerosol model：气溶胶模型，根据实际情况选择城市或农村；

• Multispectral Settings：多波段设置

K-T反演选择默认模式：Defaults->Over-Land Retrieval standard（600:2100），自动选择对应的波段，其他参数默认；

5.去云(云掩膜)

云掩膜：提取出被云覆盖的区域不参与运算

1）利用UGCS自带的云掩膜文件处理

UGCS云掩膜文件下载地址：https://glovis.usgs.gov/app?fullscreen=1

操作：打开下载的掩模图像，Toolboxs—Feature Extraction—Segmentation Image，选中云掩膜

文件，将背景值设为0

保存掩膜文件：File —Save As—Save As，选择多光谱数据—Mask，选择云检测结果作为掩膜文件，选择相应的输出的格式和路径，点击OK

得到的图像有云区域的像素值会变成nodata，当对这个图像进行快速大气校正时候，有云区域的像素将不参与计算，提高大气校正或分类精度。

2）云自动检测工具生成云掩膜（只针对Landsat8传感器）

File—Open，打开多波段数据、热红外数据和卷云波段；

辐射定标：Toolbox–Radiometric Correction(辐射校正)–Radiometric Calibration(辐射定标)

• 多光谱数据定标为大气表观反射率（TOA）

• 热红外数据（波段10/11）定标为亮度温度

• 卷云波段（波段9）定标为大气表观反射率TOA

Toolbox—Feature Extraction—Calculate Cloud Mask Using Fmask Algorithm，

Toolboxs—Feature Extraction—Segmentation Image，选中云掩膜文件，将背景值设为0；

保存掩膜文件：File —Save As—Save As，选择多光谱数据—Mask，选择云检测结果作为掩膜文件，选择相应的输出的格式和路径，点击OK；

得到的图像有云区域的像素值会变成nodata，当对这个图像进行快速大气校正时候，有云区域的像素将不参与计算，提高大气校正或分类精度。

6.图像融合

目的：将低分辨率的多光谱影像与高分辨率的单波段影像重采样生成一幅高分辨率多光谱影像，使得处理后的影像既有较高的空间分辨率，又具有多光谱特征。

操作：从File—Open中打开要融合的数据，此处以MTL.txt和band8.TIF为例。

Toolbox—Image Sharpening—Gram-Schmidt Pan Sharpening，先选择分辨率低的多光谱图像，点击OK；再选择分辨率高的单波段图像，点击OK

在参数面板中根据数据选择传感器类型，此处为：Landsat8 oli；根据需要选择重采样方法，常用重采样方法如下；最后设置输出路径和文件名。

重采样方法：

方法	概念	特点
IHS 变换	指将标准RGB图像有效地分离为代表变间信息的明度 (I)和代表波谱信息的色调 (H)、饱和度 (S)	纹理改善，空间保持较好。光谱信息损失较大，受波段限制大
Brovey	一种来自于不同传感器的数据进行融合的方法，通过归一化后的3个波段多光谱影像与高分辨率影像乘积来获得最终融合影像。	变换光谱信息保持较好，但受波段限制
乘积运算	一种直接将不同空间分辨率的影像上对应像素灰度值进行乘积运算，从而获得新的影像对应像素灰度值的融合算法。	对大的地貌类型效果好，同时可用于多光谱与高光谱的融合
PCA 变换	又称霍特林变换或K-L变换,是一种基于信息量的正交线性变换，该变换主要是采用线性投影的方法将数据投影到新的坐标空间中，从而使得新的成分按信息量分布,第一主成分包含的信息量最大，变换后各主成分分量彼此不相关，且随着主成分编号的增加该分量包含的信息量减小	无波段限制，光谱保持好，第一主成分信息高度集中，色调发生较大变化。
GS	通过统计分析方法对参与融合的各波段进行最佳匹配	改进了 PCA 中信息过分集中的问题，不受波段限制，较好的保持空间纹理信息，尤其能高保真保持光谱特征。专为最新高空间分辨率影像设计，能较好保持影像的纹理和光谱信息。

4.图像裁剪

7.图像裁剪

PS：影响裁剪之后会导致头文件丢失，因此在进行完辐射定标和大气校正之后再进行影响裁剪。

1）规则裁剪

File—New—Region of Interest，在ROI工具中选择几何工具，划定范围进行裁剪

在ROI工具箱中选择File—Save As，设置输出路径并保存；

在Toolbox中找到Regions of Interest—Subset data from ROIs，双击打开

选择要裁剪的影像—选择ROI，设置输出路径—点击OK

2）不规则裁剪——矢量数据裁剪遥感影像

打开矢量数据，打开要裁剪的遥感影像。

在Toolbox中找到Regions of Interest—Subset data from ROIs，双击打开

选择要裁剪的影像—选择ROI，设置输出路径—点击OK

8.图像镶嵌

1）打开要镶嵌的遥感影像

打开工具箱，找到Mosaicking—Seamless Mosaic，双击打开Seamless Mosaic控制面板

点击Seamless Mosaic面板中的＋号，点击select all选中需要镶嵌的数据，点击OK；

2）匀色处理—直方图匹配：

在Seamless Mosaic面板 Color Correction模块中，选择Histogram Matching，

Overlap Area Only：重叠区直方图匹配

Entire Scene：整景影像直方图匹配

3）设置透明值

在Seamless Mosaic面板的Main模块的Data Ignore Value 列表中设置透明值，当重叠区有背景值时候，可设置此值。 选择右上角的 Show Preview，可以预览镶嵌效果

4）Color Matching Action

在Seamless Mosaic面板的Main模块中,鼠标在 Color Matching Action 上右键单击，设置参考（Reference）和校正（Adjust），根据预览效果确定参考图像。设置完成后点击Finish。

5）接边线与羽化

接边线包括自动绘制和手动绘制。

在Seamless Mosaic面板的Seamlines下拉列表中选择Auto Generate Seamlines，，自动绘制接边线，自动裁剪掉 TM 边缘“锯齿”；

6）输出结果

在Seamless Mosaic面板的Export 模块中，设置参数后点击Finish。

8.波段合成

Toolbox—Raster Management—Layer Stacking，打开图层叠加面板

点击Reorder Files依据波段对图片进行排序，点击OK

三、水体提取

1.查看水体光谱曲线

看是否相对符合水体反射光谱曲线

2.计算水体指数

归一化水体提取指数              NDWI=（绿光-近红外）/（绿光+近红外）；

改进的归一化差异水体指数   MNDWI=（绿光-中红外）/（绿光+中红外）；

归一化植被差值指数             NDVI= （近红外-红光）/（近红外+红光）；

归一化植被差值指数取反      INDVI=（红光-近红外）/（红光+近红外）；

增强型植被指数                    EVI=2.5×（近红外-红光）/（近红外+6×红光-7.5×蓝光+1）

操作：Toolbox— Band Algebra—Band Math，双击打开Band Math，输入计算公式（数字应为浮点型）；

得到水体指数图：

右击layer manager中生成的数据，打开Quick Stats,发现结果范围在-1~1之间，可以进行提取；

计算阈值：绘制多个ROI来确定水体提取阈值

操作：右击数据，选自新建ROI；绘制完多边形后，右击ROI选择计算，查看阈值；

划定多个区域，直至选择出合适的阈值。

右击New Raster Color Slice，打开Edit Raster Color Slice表，编辑阈值

3.导出数据

①导出为SHP格式

点击保存—保存为Export as Class Imange，导出class image；

将水体导出成矢量：工具箱Classification——post classification——classification to vector，选择刚刚生成的二值化结果，得到水域提取的EVF数据

EVF数据转SHP：

工具栏Vector—Classic EVF to Shapefile，导入EVF数据，导出为SHP

②导出为TIFF格式

点击保存—保存为Export as Class Imange，导出class image；

右击生成的color slice文件，导出TIFF。