基于ENVI的面向坡向的地形校正方法实现
分享
-
2016-05-31
地形校正的目的
地形校正的目的主要是补偿由于不规则的地形起伏而造成的地物亮度的变化。由于这种变化会导致相似或同种植被的反射率不一致而影响遥感影像的分类精度,因此,精确的地形校正不仅能提高影像分类的精度,而且还是遥感应用的前提。
实验方法：
1、数据准备
试验中使用的影像数据是云南省的Landsat5 TM 影像, 影像获取时间为2006年01月25日, 影像中心位置位于东经99.59E,北纬28.36N,太阳高度角和方位角分别为36.15和146.46。文中选取了一块2014 pixels ×934 pixels的试验区域,该地区的地物类型主要是山林。
2、数据预处理
2.1 数据定标处
图2.1Landsat TM元数据以及DEM、NUM数据列表
如上图所示，landsat TM数据总共包含了7个波段的数据，其中6波段为热红外波段，其他6个波段为可见光波段。在软件中逐一打开7个波段，再进行波段组合给用户带来了很大的不便。利用envi主菜单 > open image file > L5131041_04120060125_MTL.txt打开文件，该文件数据可自动进行波段组合以及分类，显示结果如下图2.2：
图2.2 波段自动组合分类
波段数据读入以后进行数据的辐射定标处理，在envi下有一个专门针对Landsat数据的定标处理模块，具体的操作步骤如下：
ENVI主菜单 > Basic Tools > Preprocessing > Calibration Utilities > Landsat Calibration。
图2.3 Landsat定标模块自动读取参数
2.2 数据裁剪
数据裁剪的方式有很多，在本次试验中主要采用数据的规则采用——ROI rectangle裁剪。打开影像数据，在其image窗口中点击：Overlay > Region of Interest，弹出ROI TOOL对话框：
图2.4 ROI裁剪对话框
感兴趣区划好之后进行影像数据的裁剪：ENVI主菜单 > Basic Tools > Subset Data via ROIs。
对于DEM数据而言，该数据的投影方式以及空间分辨率可能会与上述影像数据有差，为了能将二者的数据进行很好的吻合，在进行ROI裁剪的时候需要对其ROI进行一个Map的转换，具体转换方式如下：
图2.5 ROI在不同影像间的地图转换
在此之后选择DEM数据，并对其进行裁剪，得到本实验所需要数据。
2.3 坡度、坡向数据获取
在进行地形校正的过程中，坡度、坡向对于其反应地表真实的影响很大。在获得DEM数据以后，通过ENVI主菜单下 > Topographic > Topographic Modeling进行坡度、坡向数据的提取。
图2.6 坡度、坡向提取
在地形校正过程中影响最大的坡向为南坡，需要通过BandMath工具进行南坡数据(157.5,202.5)的提取以及对应的坡度数据提取，为后面的地形校正提供数据。具体的实现方法、步骤如下：
ENVI主菜单 > Basic Tools > Band Math > (b1 ge 157.5 and b1 le 202.5)*b1+(b1 lt 157.5 and b1 gt 202.5)*0，Add to List ，点击OK，配置b1对应的波段数据Aspect文件，即可完成南坡信息提取。
图2.7 南坡信息提取条件
关于南坡影像信息对应的坡度影像，TM影像区域信息的提取方法依然采用Band Math进行提取，提取公式如上图2.7中的(b1 ne 0)*b2+(b1 eq 0)*0，其中b1代表南坡影像信息，b2代表TM影像以及坡度影像信息。重复操作两次即可提取相对于的南坡区域的TM影像以及坡度影像。
3、地形校正——坡度匹配技术
C校正模型的基本思想是:对于任意波段影像的像素DN值和其对应的太阳入射角余弦值都遵循线性关系。理想情况下,当太阳入射角为零或小于零时,表明该点缺乏太阳光照,则该点的DN值应该为零,该拟合直线应通过原点。然而,实际情况是,由于大气散射和地表相邻点反射光折射的缘故,使像素DN值和太阳入射角α成一定的余弦关系。本文采用的模型是在二阶校正模型的基础上改进建立的，需经过二个阶段（二次校正）才能得到校正结果。
按照公式（1）计算太阳有效入射角余弦值
cos i= cosz*cosS+sinz*sinS*cos(ψx-ψn) （1）
式中，z为太阳天顶角(即，90-太阳高度角)，ψx为太阳方位角，S为坡度，ψn为坡向。
在ENVI主菜单 > Basic Tools > Band Math 中输入：
cos((90-36.1530740)*0.01745)*cos(b2*0.01745)+sin((90-36.1530740) *0.01745)*sin(b2*0.01745)*cos((146.4600750-b4) *0.01745)
其中，b2 代表坡度；b4代表坡向。
第一阶段（第一次校正），按公式（2）进行。
LH = LT +（LTmax-LTmin）*（（µ-X）/µ） （2）
式中，LH为第一次校正结果，LT为原始影像，LTmax为原始影像最大值，LT民为原始影像最小值，µ为拉伸为0-255的阳坡（南坡）太阳有效入射角余弦值的平均值，X为拉伸为0-255的太阳有效入射角余弦值。
在ENVI主菜单 > Basic Tools > Band Math 中输入：
b1+(max(b1)-min(b1))*((mean(b2)-b3)/mean(b2))
其中，b1 代表Landsat 5TM原始数据；b2代表拉伸为0-255的阳坡（南坡）太阳有效入射角余弦值的平均值；b3代表拉伸为0-255的太阳有效入射角余弦值。
第二阶段（第二次校正），根据公式（3）求算模型修正系数。
Cλ=（S′λ-Nλ）/（N′λ-Nλ） （3）
式中：Cλ为图像各波段的模型修正系数，S′λ为第一次校正结果阳坡（南坡）平均值，N′λ为第一次校正结果阴坡（北坡）平均值，Nλ为原始图像阴坡（北坡）平均值。
在ENVI主菜单 > Basic Tools > Band Math 中输入：
(mean(b1)-mean(b2))/(mean(b3)-mean(b2))
其中，b1 代表第一次校正结果阳坡（南坡）平均值；b2代表原始图像阴坡（北坡）平均值；b3代表第一次校正结果阴坡（北坡）平均值。
把计算出的各波段模型修正系数带入公式（2）建立公式（4）进行第二次校正。
LH = LT +（LTmax-LTmin）*（（µ-X）/µ）* Cλ （4）
在ENVI主菜单 > Basic Tools > Band Math 中输入：
b1+(max(b2)-min(b2))*((225.43-b3)/225.43)* (mean(b4)-mean(b5))/(mean(b6)-mean(b5))
其中，其中，b1 代表Landsat 5TM原始数据；b2代表拉伸为0-255的阳坡（南坡）太阳有效入射角余弦值的平均值；b3代表拉伸为0-255的太阳有效入射角余弦值。b4 代表第一次校正结果阳坡（南坡）平均值；b5代表原始图像阴坡（北坡）平均值，b6第一次校正结果阴坡（北坡）平均值。
4、结果展示
图 4.1 原始影像、第一次校正结果以及坡向校正结果对比图
图 4.2 原始影像、第一次校正结果以及坡向校正结果对比图
根据上述两幅图的对比可以看出，在经过坡向地形校正以后其校正影像显示结果很显著。需要注意的是：利用该种方法进行地形校正的时候一定要保证dem数据的准确性，不然在后面进行坡度坡向信息提取的时候容易出现与实际情况不符的状况，最终影像到纠正结果。
结果分析
采用该种方法很好的进行了地形校正，是不是最终的结果能将山体阴影下的植被给予很好的显示，需要对其校正前后获取的NDVI值进行对比分析，从而确定该种方法的精度。
图5.1 校正前后的NDVI值对比
图5.2 山体阴影处经过校正后的NDVI值的变化
上述两幅图是从山体阴影上分析校正前后NDVI值的变化情况，从这个DN值的变化上可以看出，经过校正后的影像获取的NDVI值更为接近实际地表的植被覆盖情况。对于裸露的地表，其NDVI值的变化是不是很大，下面从裸土着手分析，坡向地形校正对于不存在阴影的地区的校正结果如何。
图 5.3裸露地表NDVI值的变化
从上图5.3的结果对比分析可以看出，裸露地表的NDVI值的变化在0.03范围，变化不是很大，这与开始设想的坡向地形校正方法主要是针对山体阴影区域的校正，对于平坦区域的作用影像还是比较小的。
综合上述的结果对比分析以及影像显示效果分析得出，该方法能很好的帮助我们进行山体阴影区域的植被生长情况的恢复，为以后的研究工作提供了更为便捷的校正方法。
致谢： 本文所用的数据和方法由西南林业大学林业3S技术工程中心何超老师提供。在此感谢何老师的支持与帮助！
文章来源：http://blog.sina.com.cn/s/blog_764b1e9d0100vax9.html
地形校正的目的主要是补偿由于不规则的地形起伏而造成的地物亮度的变化。由于这种变化会导致相似或同种植被的反射率不一致而影响遥感影像的分类精度,因此,精确的地形校正不仅能提高影像分类的精度,而且还是遥感应用的前提。
实验方法：
1、数据准备
试验中使用的影像数据是云南省的Landsat5 TM 影像, 影像获取时间为2006年01月25日, 影像中心位置位于东经99.59E,北纬28.36N,太阳高度角和方位角分别为36.15和146.46。文中选取了一块2014 pixels ×934 pixels的试验区域,该地区的地物类型主要是山林。
2、数据预处理
2.1 数据定标处
图2.1Landsat TM元数据以及DEM、NUM数据列表
如上图所示，landsat TM数据总共包含了7个波段的数据，其中6波段为热红外波段，其他6个波段为可见光波段。在软件中逐一打开7个波段，再进行波段组合给用户带来了很大的不便。利用envi主菜单 > open image file > L5131041_04120060125_MTL.txt打开文件，该文件数据可自动进行波段组合以及分类，显示结果如下图2.2：
图2.2 波段自动组合分类
波段数据读入以后进行数据的辐射定标处理，在envi下有一个专门针对Landsat数据的定标处理模块，具体的操作步骤如下：
ENVI主菜单 > Basic Tools > Preprocessing > Calibration Utilities > Landsat Calibration。
图2.3 Landsat定标模块自动读取参数
2.2 数据裁剪
数据裁剪的方式有很多，在本次试验中主要采用数据的规则采用——ROI rectangle裁剪。打开影像数据，在其image窗口中点击：Overlay > Region of Interest，弹出ROI TOOL对话框：
图2.4 ROI裁剪对话框
感兴趣区划好之后进行影像数据的裁剪：ENVI主菜单 > Basic Tools > Subset Data via ROIs。
对于DEM数据而言，该数据的投影方式以及空间分辨率可能会与上述影像数据有差，为了能将二者的数据进行很好的吻合，在进行ROI裁剪的时候需要对其ROI进行一个Map的转换，具体转换方式如下：
图2.5 ROI在不同影像间的地图转换
在此之后选择DEM数据，并对其进行裁剪，得到本实验所需要数据。
2.3 坡度、坡向数据获取
在进行地形校正的过程中，坡度、坡向对于其反应地表真实的影响很大。在获得DEM数据以后，通过ENVI主菜单下 > Topographic > Topographic Modeling进行坡度、坡向数据的提取。
图2.6 坡度、坡向提取
在地形校正过程中影响最大的坡向为南坡，需要通过BandMath工具进行南坡数据(157.5,202.5)的提取以及对应的坡度数据提取，为后面的地形校正提供数据。具体的实现方法、步骤如下：
ENVI主菜单 > Basic Tools > Band Math > (b1 ge 157.5 and b1 le 202.5)*b1+(b1 lt 157.5 and b1 gt 202.5)*0，Add to List ，点击OK，配置b1对应的波段数据Aspect文件，即可完成南坡信息提取。
图2.7 南坡信息提取条件
关于南坡影像信息对应的坡度影像，TM影像区域信息的提取方法依然采用Band Math进行提取，提取公式如上图2.7中的(b1 ne 0)*b2+(b1 eq 0)*0，其中b1代表南坡影像信息，b2代表TM影像以及坡度影像信息。重复操作两次即可提取相对于的南坡区域的TM影像以及坡度影像。
3、地形校正——坡度匹配技术
C校正模型的基本思想是:对于任意波段影像的像素DN值和其对应的太阳入射角余弦值都遵循线性关系。理想情况下,当太阳入射角为零或小于零时,表明该点缺乏太阳光照,则该点的DN值应该为零,该拟合直线应通过原点。然而,实际情况是,由于大气散射和地表相邻点反射光折射的缘故,使像素DN值和太阳入射角α成一定的余弦关系。本文采用的模型是在二阶校正模型的基础上改进建立的，需经过二个阶段（二次校正）才能得到校正结果。
按照公式（1）计算太阳有效入射角余弦值
cos i= cosz*cosS+sinz*sinS*cos(ψx-ψn) （1）
式中，z为太阳天顶角(即，90-太阳高度角)，ψx为太阳方位角，S为坡度，ψn为坡向。
在ENVI主菜单 > Basic Tools > Band Math 中输入：
cos((90-36.1530740)*0.01745)*cos(b2*0.01745)+sin((90-36.1530740) *0.01745)*sin(b2*0.01745)*cos((146.4600750-b4) *0.01745)
其中，b2 代表坡度；b4代表坡向。
第一阶段（第一次校正），按公式（2）进行。
LH = LT +（LTmax-LTmin）*（（µ-X）/µ） （2）
式中，LH为第一次校正结果，LT为原始影像，LTmax为原始影像最大值，LT民为原始影像最小值，µ为拉伸为0-255的阳坡（南坡）太阳有效入射角余弦值的平均值，X为拉伸为0-255的太阳有效入射角余弦值。
在ENVI主菜单 > Basic Tools > Band Math 中输入：
b1+(max(b1)-min(b1))*((mean(b2)-b3)/mean(b2))
其中，b1 代表Landsat 5TM原始数据；b2代表拉伸为0-255的阳坡（南坡）太阳有效入射角余弦值的平均值；b3代表拉伸为0-255的太阳有效入射角余弦值。
第二阶段（第二次校正），根据公式（3）求算模型修正系数。
Cλ=（S′λ-Nλ）/（N′λ-Nλ） （3）
式中：Cλ为图像各波段的模型修正系数，S′λ为第一次校正结果阳坡（南坡）平均值，N′λ为第一次校正结果阴坡（北坡）平均值，Nλ为原始图像阴坡（北坡）平均值。
在ENVI主菜单 > Basic Tools > Band Math 中输入：
(mean(b1)-mean(b2))/(mean(b3)-mean(b2))
其中，b1 代表第一次校正结果阳坡（南坡）平均值；b2代表原始图像阴坡（北坡）平均值；b3代表第一次校正结果阴坡（北坡）平均值。
把计算出的各波段模型修正系数带入公式（2）建立公式（4）进行第二次校正。
LH = LT +（LTmax-LTmin）*（（µ-X）/µ）* Cλ （4）
在ENVI主菜单 > Basic Tools > Band Math 中输入：
b1+(max(b2)-min(b2))*((225.43-b3)/225.43)* (mean(b4)-mean(b5))/(mean(b6)-mean(b5))
其中，其中，b1 代表Landsat 5TM原始数据；b2代表拉伸为0-255的阳坡（南坡）太阳有效入射角余弦值的平均值；b3代表拉伸为0-255的太阳有效入射角余弦值。b4 代表第一次校正结果阳坡（南坡）平均值；b5代表原始图像阴坡（北坡）平均值，b6第一次校正结果阴坡（北坡）平均值。
4、结果展示
图 4.1 原始影像、第一次校正结果以及坡向校正结果对比图
图 4.2 原始影像、第一次校正结果以及坡向校正结果对比图
根据上述两幅图的对比可以看出，在经过坡向地形校正以后其校正影像显示结果很显著。需要注意的是：利用该种方法进行地形校正的时候一定要保证dem数据的准确性，不然在后面进行坡度坡向信息提取的时候容易出现与实际情况不符的状况，最终影像到纠正结果。
结果分析
采用该种方法很好的进行了地形校正，是不是最终的结果能将山体阴影下的植被给予很好的显示，需要对其校正前后获取的NDVI值进行对比分析，从而确定该种方法的精度。
图5.1 校正前后的NDVI值对比
图5.2 山体阴影处经过校正后的NDVI值的变化
上述两幅图是从山体阴影上分析校正前后NDVI值的变化情况，从这个DN值的变化上可以看出，经过校正后的影像获取的NDVI值更为接近实际地表的植被覆盖情况。对于裸露的地表，其NDVI值的变化是不是很大，下面从裸土着手分析，坡向地形校正对于不存在阴影的地区的校正结果如何。
图 5.3裸露地表NDVI值的变化
从上图5.3的结果对比分析可以看出，裸露地表的NDVI值的变化在0.03范围，变化不是很大，这与开始设想的坡向地形校正方法主要是针对山体阴影区域的校正，对于平坦区域的作用影像还是比较小的。
综合上述的结果对比分析以及影像显示效果分析得出，该方法能很好的帮助我们进行山体阴影区域的植被生长情况的恢复，为以后的研究工作提供了更为便捷的校正方法。
致谢： 本文所用的数据和方法由西南林业大学林业3S技术工程中心何超老师提供。在此感谢何老师的支持与帮助！
文章来源：http://blog.sina.com.cn/s/blog_764b1e9d0100vax9.html
