ENVI支持下的遥感荒漠化地信息提取
分享
-
2016-06-18
背景
土地荒漠化是指包括气候变异和人类活动在内的种种因素造成的干旱半干旱和亚湿润干地区的土地退化。及时准确地掌握土地荒漠化发生发展情况是有效防止和治理土地荒漠化的基本前提。目前遥感技术在土地荒漠化监测中起到了不可替代的作用。使用遥感影像数据可以提取土地荒漠化信息，通过遥感影像所表现的不同信息，可以判断土地荒漠化的发生与否以及发展程度等。在进行土地荒漠化信息提取时，常用的方法有人工目视解译方法、监督分类方法、非监督分类方法、决策树分层分类方法、神经网络自动提取方法等。在实际应用中，通常选择其中的一种或结合几种方法进行分类提取。
目前，一种比较新的方法是通过构造“植被指数(NDVI)——反照率(Albedo)特征空间”来进行荒漠化信息遥感提取。荒漠化过程及其地表特性的变化能在 Albedo-NDVI特征空间中得到明显直观的反映。在Albedo-NDVI特征空间中，可以利用植被指数和地表反照率的组合信息，通过选择反映荒漠化程度的合理指数，就可以将不同荒漠化土地有效地加以区分，从而实现荒漠化时空分布与动态变化的定量监测与研究。而这个问题的合理解决，实际上就是如何根据需要采用一定的综合指标来划分Albedo-NDVI特征空间。
根据Verstraete and Pinty的研究结论，如果在代表荒漠化变化趋势的垂直方向上划分Albedo-NDVI特征空间，可以将不同的荒漠化土地有效地区分开来。而垂线方向在 Albedo-NDVI特征空间的位置可以用特征空间中简单的二元线性多项式加以表达，即DDI=k*NDVI - Albedo，其中: DDI可称为荒漠化分级指数，k由特征空间中拟合的曲线斜率确定。
本文主要介绍在ENVI下实现Albedo-NDVI特征空间遥感荒漠化信息提取的操作流程。
2．处理流程介绍
流程说明：
（1）数据获取
（2）数据预处理
包括数据定标处理、大气校正、几何配准、研究区域的裁剪。
（3）信息提取
根据前人研究的公式计算NDVI和Albedo。然后将结果进行归一化处理，保证数据的一致性。
（4）计算NDVI和Albedo的定量关系
这一步利用excel工具找到NDVI和Albedo数据间的量化关系。得到表达式Albedo=a*NDVI+b中a的值。
（5）荒漠化差值指数的计算
利用bandmath工具实现表达式：DDI=(-1/a)*NDVI-Albedo，得到荒漠化差值植被指数。
（6）荒漠化分级信息的提取
根据荒漠化差值植被指数就能进行荒漠化分级信息提取。有两种方法：一是通过设置分级阈值进行分级；另一是通过利用“自然间断点分级法”将DDI值进行分级。
（7）制图输出
将结果图输出。
3．详细处理过程
以下有部分基础内容没有详细描述，在阅读这些内容时，假设您已经具备了基本的ENVI操作知识。
3.1 数据预处理
第一步：数据读取和定标
主菜单->File->Open External File->Landsat->Geo TIFF with Metadata，打开TM数据L51124038_03820020902_MTL.txt
第二步：数据定标
主菜单->Basic Tools->Preprocesssing->Calibration Utilities->Landsat Calibration，选择多波段文件，弹出ENVI Landsat Calibration对话框，自动读取元数据中的信息并加载（如图），定标类型选择Radiance，指定保存路径，点击OK。
第三步：几何校正和裁剪
这一步是对环境数据进行几何校正，使其具有精确的地理信息，基准影像可以是地形图，也可以是已经过校正的其他中高分辨率影像，本实验中用到的是已经过校正的TM数据作为基准影像。由于整景数据范围非常大，而我们的工程区只是其中一小部分，所以在进行几何配准之前，将研究区域裁剪出来。
第四步：大气校正
（一）FLAASH对图像文件有以下几个要求：
1) 数据是经过定标后的辐射亮度（辐射率）数据，单位是：（μW）/（cm2*nm*sr）。
2) 数据带有中心波长（wavelenth）值，如果是高光谱还必须有波段宽度（FWHM）,这两个参数都可以通过编辑头文件信息输入（Edit Header）。
3) 数据类型
支持四种数据类型：浮点型（floating）、长整型(long integer )、整型（integer）和无符号整型 (unsignedint)。数据存储类型： ENVI标准栅格格式文件，且是BIP或者BIL。
4) 波谱范围：0.4－2500μm。
本次用的TM数据经过以上处理，如已经定标为W*m^(-2)*sr^(-1)*um^(-1)单位、浮点型的辐射率数据，有中心波长信息，下面将BSQ格式转成BIL格式。
选择主菜单Basic Tools->Convert Data（BSQ、BIL、BIP），选择已经经过配准的TM-Cal-sub-jz.img，在Convert File Parameters中，Output Interleave选择BIL，选择Convert In Place：yes，单击OK，
（如下图）。
（二）：设置参数进行FLAASH大气校正
（1） 主菜单Spectral->FLAASH打开FLAASH大气校正模块；
（2） 点击Input Radiance Image，选择BIL格式的TM-Cal-sub-jz.img，在Radiance Scale Factors面板中选择Use single scale factor for all bands，由于定标的辐射量数据与FLAASH的辐射亮度的单位相差10倍，所以在此Single scale factor选择默认：10，单击OK；
（3） 设置输出文件及路径设置
（4） 传感器基本信息设置：
成像中心点经纬度、传感器高度、成像区域平均高度、成像时间，设置，这些都可以从数据头文件中读取L71124038_03820020902_MTL.txt。
（5） 大气模型，选择MLS，气溶胶模型，选择Rural，气溶胶反演方法选择2-Band（K-T），能见度给40km。
（6） 单击Advanced Settings，在高级设置中，Tile Size 默认的是Cash size 的大小，手动改为100Mb，单击OK；
（7） 设置好后，在大气校正模块面板中，单击Apply。
（8） 大气校正完成后，检查大气校正的结果，分别加载校正前后的图像，将两幅影像进行地理链接，移动到植被区域（植被的波谱曲线比较特殊），在影像上右键，选择Z Profile（Spectrum）打开光谱曲线窗口，显示两幅图像同一位置的光谱曲线图
3.2 信息提取
第一步：提取NDVI
主菜单->Transform->NDVI，弹出NDVI Calculation Parameters对话框，自动标定波段。保存文件，输出，点击OK。
第二步：提取Albedo
地表反照率反演：利用 L i a n g 建立的 L a n d s a t — T M数据的反演模型，估算了研究区地表反照率。
（1） 主菜单->Basic Tools->Bandmath->弹出Band Math对话框，键入表达式：0.356*b1+0.13*b3+0.373*b4+0.085*b5+0.072*b7-0.0018，点击Add to List，点击OK。
（2） 在弹出的Variables to Band Pairings对话框中分别为B1、B3…、B7指定相应的波段（经过大气校正后的数据）。
第三步：归一化处理
采用归一化公式进行NDVI和albedo数据的归一化处理。归一化公式如下。
（1）计算NDVI和albedo数据的最大最小值。通过波段列表（Available Band List）中选择NDVI和albedo文件，右键该文件，点击Quick Stats…进行统计后，弹出Statistics Results对话框，可以获取NDVI或者albedo的最大最小值。
（2）归一化计算。主菜单->Basic Tools->Bandmath->弹出Band Math对话框，键入表达式：(b1+0.429861)/(0.742756+0.429861)，点击Add to List，点击OK。在弹出的Variables to Band Pairings对话框中分别为B1指定相应的波段（NDVI数据）。
（3）按照同样的方法，归一化Albedo数据。
3.3 计算NDVI与albedo的定量关系
经专家学者研究NDVI与albedo之间存在着一种负相关的线性关系。不同沙漠化土地类型对应的植被指数 ( N D V I ) 和地表反照率 ( A l b e d o ) 具有显著的线性负相关性，类似于：Albedo=65.324-0.453*NDVI这种形式，通用表达式Albedo=a*NDVI+b。说明，随着荒漠化程度的增加，植被指数 ( N D V I ) 逐渐减少，而地表反照率则逐渐增加。在 A l b e d o—NDVI 特征空间中，荒漠化过程得到了明显的反映。
为了找到两者之间的定量关系，需要分别找出NDVI和albedo对应的两组数据，利用这两组数据进行回归拟合出一个关系式。
第一步：选择相应点
（1）在NDVI或者Albedo的图像窗口中，右键->选择ROI Tools，弹出ROI Tools对话框在ROI_Type中选择Point。然后点击Image，在image窗口中选点。
（2）选好点后，将点导出。
在ROI TOOL中，选择File->Output ROIs to ASCII。选择NDVI的图像，在Output ROIs to ASCII Parameters面板中，选择ROI点，单击Edit Output ASCII Form，在输出内容设置面板中（如下图），选择ID、经纬度（Geo Location）、和波段像元值（Band Values）。点击OK。指定输出路径和名称，点击OK，将对应的NDVI点值输出。同样的方法前面选择的ROI点对应的Albedo的点值输出为Albedo.txt文件。
第二步：计算定量关系
在EXCEL软件中进行线性拟合两者的定量关系。
有了相同位置的NDVI值和Albedo值，在Excel中选中ndvi值与albedo值，绘制散点图。
在散点图上选中散点，单击右键->添加趋势线，打开设置趋势线格式面板，勾选线性，显示公式，显示R平方值。点击“关闭”按钮，线性回归方程和R方值在散点图上显示。
添加趋势线
线性反演模型
得到反演模型为：y=-5.062x+5.182，R2=0.841。
3.4荒漠化差值指数的计算
通过上一步的处理得到了参数a，根据公式a*k=-1，可以计算出k。将k值带入荒漠化差值指数表达式
中可以计算DDI。表达为：DDI=1/5.062*NDVI-Albedo，使用Basic Tool->Bandmath，在Enter an expression下面输入表达式：1/5.062*b1 –b2，单击Add to List，单击OK，在Variables to Bands Pairings面板中，选择b1为NDVI的图像，b2为Albedo的图像，设置输出路径和文件名，单击OK，计算得到荒漠化差值指数的反演图。
3.5荒漠化分级信息提取
方法一：
实地考察，根据相关标准，将该区域的荒漠化程度分级，即非荒漠化、轻度荒漠化、中度荒漠化、重度荒漠化和极重度荒漠化。找出不同荒漠化级别与对应的荒漠化差值指数图上的临界点。然后利用Density Slice工具进行分级显示。
在Display中显示荒漠化差值指数，是一个灰度的单波段图像。
(1) 选择Tools->Color Mapping->Density Slice，单击Clear Range按钮清除默认区间。
(2) 选择Opions->Add New Ranges，增加四个区间
(3) 单击Apply。
(4) 选择File->Output Range to Class Image，可以将反演结果输出。
方法二：
利用自然间断点分级法( natura l break)将DDI值进行分级，这种分类方法是基于统计学的Jenk最优化法得出的分界点，能够使各级的内部方差之和最小。有不少人采用这种方法。
将DDI图像数据加载到ArcMAP中，打开工具箱ArcToolBox->空间分析工具（Spatial Analyst）->重分类（Reclassify）->重分类。弹出重分类对话框如图。选中数据、点击“分类”按钮，弹出分类对话框，选择“自然间断点分级法”、设设置“类别”数，点击“确定”。
文章来源：http://blog.sina.com.cn/s/blog_764b1e9d0100ywdl.html
土地荒漠化是指包括气候变异和人类活动在内的种种因素造成的干旱半干旱和亚湿润干地区的土地退化。及时准确地掌握土地荒漠化发生发展情况是有效防止和治理土地荒漠化的基本前提。目前遥感技术在土地荒漠化监测中起到了不可替代的作用。使用遥感影像数据可以提取土地荒漠化信息，通过遥感影像所表现的不同信息，可以判断土地荒漠化的发生与否以及发展程度等。在进行土地荒漠化信息提取时，常用的方法有人工目视解译方法、监督分类方法、非监督分类方法、决策树分层分类方法、神经网络自动提取方法等。在实际应用中，通常选择其中的一种或结合几种方法进行分类提取。
目前，一种比较新的方法是通过构造“植被指数(NDVI)——反照率(Albedo)特征空间”来进行荒漠化信息遥感提取。荒漠化过程及其地表特性的变化能在 Albedo-NDVI特征空间中得到明显直观的反映。在Albedo-NDVI特征空间中，可以利用植被指数和地表反照率的组合信息，通过选择反映荒漠化程度的合理指数，就可以将不同荒漠化土地有效地加以区分，从而实现荒漠化时空分布与动态变化的定量监测与研究。而这个问题的合理解决，实际上就是如何根据需要采用一定的综合指标来划分Albedo-NDVI特征空间。
根据Verstraete and Pinty的研究结论，如果在代表荒漠化变化趋势的垂直方向上划分Albedo-NDVI特征空间，可以将不同的荒漠化土地有效地区分开来。而垂线方向在 Albedo-NDVI特征空间的位置可以用特征空间中简单的二元线性多项式加以表达，即DDI=k*NDVI - Albedo，其中: DDI可称为荒漠化分级指数，k由特征空间中拟合的曲线斜率确定。
本文主要介绍在ENVI下实现Albedo-NDVI特征空间遥感荒漠化信息提取的操作流程。
2．处理流程介绍
流程说明：
（1）数据获取
（2）数据预处理
包括数据定标处理、大气校正、几何配准、研究区域的裁剪。
（3）信息提取
根据前人研究的公式计算NDVI和Albedo。然后将结果进行归一化处理，保证数据的一致性。
（4）计算NDVI和Albedo的定量关系
这一步利用excel工具找到NDVI和Albedo数据间的量化关系。得到表达式Albedo=a*NDVI+b中a的值。
（5）荒漠化差值指数的计算
利用bandmath工具实现表达式：DDI=(-1/a)*NDVI-Albedo，得到荒漠化差值植被指数。
（6）荒漠化分级信息的提取
根据荒漠化差值植被指数就能进行荒漠化分级信息提取。有两种方法：一是通过设置分级阈值进行分级；另一是通过利用“自然间断点分级法”将DDI值进行分级。
（7）制图输出
将结果图输出。
3．详细处理过程
以下有部分基础内容没有详细描述，在阅读这些内容时，假设您已经具备了基本的ENVI操作知识。
3.1 数据预处理
第一步：数据读取和定标
主菜单->File->Open External File->Landsat->Geo TIFF with Metadata，打开TM数据L51124038_03820020902_MTL.txt
第二步：数据定标
主菜单->Basic Tools->Preprocesssing->Calibration Utilities->Landsat Calibration，选择多波段文件，弹出ENVI Landsat Calibration对话框，自动读取元数据中的信息并加载（如图），定标类型选择Radiance，指定保存路径，点击OK。
第三步：几何校正和裁剪
这一步是对环境数据进行几何校正，使其具有精确的地理信息，基准影像可以是地形图，也可以是已经过校正的其他中高分辨率影像，本实验中用到的是已经过校正的TM数据作为基准影像。由于整景数据范围非常大，而我们的工程区只是其中一小部分，所以在进行几何配准之前，将研究区域裁剪出来。
第四步：大气校正
（一）FLAASH对图像文件有以下几个要求：
1) 数据是经过定标后的辐射亮度（辐射率）数据，单位是：（μW）/（cm2*nm*sr）。
2) 数据带有中心波长（wavelenth）值，如果是高光谱还必须有波段宽度（FWHM）,这两个参数都可以通过编辑头文件信息输入（Edit Header）。
3) 数据类型
支持四种数据类型：浮点型（floating）、长整型(long integer )、整型（integer）和无符号整型 (unsignedint)。数据存储类型： ENVI标准栅格格式文件，且是BIP或者BIL。
4) 波谱范围：0.4－2500μm。
本次用的TM数据经过以上处理，如已经定标为W*m^(-2)*sr^(-1)*um^(-1)单位、浮点型的辐射率数据，有中心波长信息，下面将BSQ格式转成BIL格式。
选择主菜单Basic Tools->Convert Data（BSQ、BIL、BIP），选择已经经过配准的TM-Cal-sub-jz.img，在Convert File Parameters中，Output Interleave选择BIL，选择Convert In Place：yes，单击OK，
（如下图）。
（二）：设置参数进行FLAASH大气校正
（1） 主菜单Spectral->FLAASH打开FLAASH大气校正模块；
（2） 点击Input Radiance Image，选择BIL格式的TM-Cal-sub-jz.img，在Radiance Scale Factors面板中选择Use single scale factor for all bands，由于定标的辐射量数据与FLAASH的辐射亮度的单位相差10倍，所以在此Single scale factor选择默认：10，单击OK；
（3） 设置输出文件及路径设置
（4） 传感器基本信息设置：
成像中心点经纬度、传感器高度、成像区域平均高度、成像时间，设置，这些都可以从数据头文件中读取L71124038_03820020902_MTL.txt。
（5） 大气模型，选择MLS，气溶胶模型，选择Rural，气溶胶反演方法选择2-Band（K-T），能见度给40km。
（6） 单击Advanced Settings，在高级设置中，Tile Size 默认的是Cash size 的大小，手动改为100Mb，单击OK；
（7） 设置好后，在大气校正模块面板中，单击Apply。
（8） 大气校正完成后，检查大气校正的结果，分别加载校正前后的图像，将两幅影像进行地理链接，移动到植被区域（植被的波谱曲线比较特殊），在影像上右键，选择Z Profile（Spectrum）打开光谱曲线窗口，显示两幅图像同一位置的光谱曲线图
3.2 信息提取
第一步：提取NDVI
主菜单->Transform->NDVI，弹出NDVI Calculation Parameters对话框，自动标定波段。保存文件，输出，点击OK。
第二步：提取Albedo
地表反照率反演：利用 L i a n g 建立的 L a n d s a t — T M数据的反演模型，估算了研究区地表反照率。
（1） 主菜单->Basic Tools->Bandmath->弹出Band Math对话框，键入表达式：0.356*b1+0.13*b3+0.373*b4+0.085*b5+0.072*b7-0.0018，点击Add to List，点击OK。
（2） 在弹出的Variables to Band Pairings对话框中分别为B1、B3…、B7指定相应的波段（经过大气校正后的数据）。
第三步：归一化处理
采用归一化公式进行NDVI和albedo数据的归一化处理。归一化公式如下。
（1）计算NDVI和albedo数据的最大最小值。通过波段列表（Available Band List）中选择NDVI和albedo文件，右键该文件，点击Quick Stats…进行统计后，弹出Statistics Results对话框，可以获取NDVI或者albedo的最大最小值。
（2）归一化计算。主菜单->Basic Tools->Bandmath->弹出Band Math对话框，键入表达式：(b1+0.429861)/(0.742756+0.429861)，点击Add to List，点击OK。在弹出的Variables to Band Pairings对话框中分别为B1指定相应的波段（NDVI数据）。
（3）按照同样的方法，归一化Albedo数据。
3.3 计算NDVI与albedo的定量关系
经专家学者研究NDVI与albedo之间存在着一种负相关的线性关系。不同沙漠化土地类型对应的植被指数 ( N D V I ) 和地表反照率 ( A l b e d o ) 具有显著的线性负相关性，类似于：Albedo=65.324-0.453*NDVI这种形式，通用表达式Albedo=a*NDVI+b。说明，随着荒漠化程度的增加，植被指数 ( N D V I ) 逐渐减少，而地表反照率则逐渐增加。在 A l b e d o—NDVI 特征空间中，荒漠化过程得到了明显的反映。
为了找到两者之间的定量关系，需要分别找出NDVI和albedo对应的两组数据，利用这两组数据进行回归拟合出一个关系式。
第一步：选择相应点
（1）在NDVI或者Albedo的图像窗口中，右键->选择ROI Tools，弹出ROI Tools对话框在ROI_Type中选择Point。然后点击Image，在image窗口中选点。
（2）选好点后，将点导出。
在ROI TOOL中，选择File->Output ROIs to ASCII。选择NDVI的图像，在Output ROIs to ASCII Parameters面板中，选择ROI点，单击Edit Output ASCII Form，在输出内容设置面板中（如下图），选择ID、经纬度（Geo Location）、和波段像元值（Band Values）。点击OK。指定输出路径和名称，点击OK，将对应的NDVI点值输出。同样的方法前面选择的ROI点对应的Albedo的点值输出为Albedo.txt文件。
第二步：计算定量关系
在EXCEL软件中进行线性拟合两者的定量关系。
有了相同位置的NDVI值和Albedo值，在Excel中选中ndvi值与albedo值，绘制散点图。
添加趋势线
线性反演模型
得到反演模型为：y=-5.062x+5.182，R2=0.841。
3.4荒漠化差值指数的计算
通过上一步的处理得到了参数a，根据公式a*k=-1，可以计算出k。将k值带入荒漠化差值指数表达式
3.5荒漠化分级信息提取
方法一：
实地考察，根据相关标准，将该区域的荒漠化程度分级，即非荒漠化、轻度荒漠化、中度荒漠化、重度荒漠化和极重度荒漠化。找出不同荒漠化级别与对应的荒漠化差值指数图上的临界点。然后利用Density Slice工具进行分级显示。
在Display中显示荒漠化差值指数，是一个灰度的单波段图像。
(1) 选择Tools->Color Mapping->Density Slice，单击Clear Range按钮清除默认区间。
(2) 选择Opions->Add New Ranges，增加四个区间
(3) 单击Apply。
(4) 选择File->Output Range to Class Image，可以将反演结果输出。
方法二：
利用自然间断点分级法( natura l break)将DDI值进行分级，这种分类方法是基于统计学的Jenk最优化法得出的分界点，能够使各级的内部方差之和最小。有不少人采用这种方法。
将DDI图像数据加载到ArcMAP中，打开工具箱ArcToolBox->空间分析工具（Spatial Analyst）->重分类（Reclassify）->重分类。弹出重分类对话框如图。选中数据、点击“分类”按钮，弹出分类对话框，选择“自然间断点分级法”、设设置“类别”数，点击“确定”。
文章来源：http://blog.sina.com.cn/s/blog_764b1e9d0100ywdl.html
