遥感影像空间增强与傅立叶变换滤波处理

1、实习序号及题目第四次遥感影像空间增强与傅立叶变换滤波处理实习人姓名薛专业班级2014级地理信息科学e-mail地址实习指导教师姓名党涛，李亚宁，任建平，韩希光实习地点榆中校区实验楼A109实习日期时间2016年 11月 18日实习目的和内容1.实习目的：理解遥感影像空间增强和傅立叶变换的概念和意义，掌握运用ENVI/IDL进行空间增强和傅立叶变换及滤波进行地形和地貌信息提取的步骤和方法。学会使用SRTM DEM数据应用的初步步骤。2.实习内容：利用下述操作对下载的DEM数据进行处理，比较变换前后影像的显示结果1. 从下载的ASTER DEM中选取512512子区数据，显示分析研究区地形地貌特征

2、2. 完成如下卷积增强运算，比较变换前后影像中地形/地貌信息上的差异：l 高通滤波l 拉普拉斯变换l 定向增强l 高斯高通滤波l Sobel边缘增强l Roberts算子增强3. 纹理运算，比较变换前后影像中地形/地貌信息上的差异：l 一阶纹理运算l 二阶纹理运算4. 周期噪声去除：利用傅立叶变换对影像tm_1.img所有7个波段进行处理，选择恰当的参数消除影像上的噪声信号原理和方法 基本概念：1. 空间卷积：又称为空间滤波（filtering），系采用空间域的处理方法，在邻域像元的参不下迚行运算处理，由此生成的影像中各个像元的值是原影像对应像元及其邻域像元灰度值的加权和。通过空间卷积运算达到

3、突出影像空间中的变化或变异性规律，或抑制噪声信号的目的。设窗口大小为mn，（i,j）是窗口中心像元，f(x,y)是影像像元值，g(i,j)是运算结果，h(x,y)是窗口模板，那么卷积计算的一般公式为：m,ng(I,j)=(f(x,y)h(x,y) y=1,x=1其中，模板即积运算的邻域像元权重系数矩阵，定制的邻域运算关系式，又称卷积核（convolution kernel）、卷积函数、滤波器（filtering）。2. 纹理分析：纹理是一种重要的视觉信息，主要体现遥感图像中局部灰度或者彩色的有规律的重复出现。它是影像中第五形状、大小、色彩的综合体现，为遥感影像的目视解译提供了有力的信息。通过一

4、阶纹理运算与二阶纹理运算可以比较影像地形、地物上的差异。3. 傅里叶变换消除噪声是在完成了傅里叶变换之后，通过观察频率域影像，选取高亮部分（图像上单位像元与邻域像元差异较大），制作成滤波模板，然后用制作好的滤波模板完成对信号的剔除，形成新的频率影像，最后通过傅里叶逆变换，完成从频率域到空间域的变换，得到去除噪声的影像。一维傅里叶变换：逆变换：二维傅里叶变换：逆变换：4. 拉普拉斯变换是工程数学中常用的一种积分变换，又名拉氏变换。拉普拉斯变换是一个线性变换，可将一个有参数实数t（t 0）的函数转换为一个参数为复数s的函数。 拉普拉斯变换是对于t=0函数值不为零的连续时间函数x(t)通过关系式 (

5、式中st为自然对数底e的指数)变换为复变量s的函数X(s)。它也是时间函数x(t)的“复频域”表示方式。5. 边缘增强能提高边缘与周围像素之间的反差。线性高通滤波器是最常用的线性锐化滤波器，中心系数为正数，其他系数为负数；梯度法是非线性锐化滤波器，在梯度图像上梯度值较大的部分就是边缘。罗伯特Roberts梯度采用交叉差分的方法，检测像素与其在上下之间或左右之间或斜方向之间的差异；Sobel 梯度对 4-邻域采用加权方法进行差分，对边缘的检测更加精确；定向检测在提取边缘时指定了方向，有目的的提取某一特定方向的边缘或线性特征。6. Sobel算子包含两组3x3的矩阵，分别为横向及纵向，将之与图像作

6、平面卷积，即可分别得出横向及纵向的亮度差分近似值。如果以A代表原始图像，Gx及Gy分别代表经横向及纵向边缘检测的图像，其公式如下:图像的每一个像素的横向及纵向梯度近似值可用以下的公式结合，来计算梯度的大小。然后可用以下公式计算梯度方向在以上例子中，如果以上的角度等于零，即代表图像该处拥有纵向边缘，左方较右方暗。数据准备与研究区概况1.数据准备：Landsat tm_1.img；ASTGTM_N35E1032.数据特征：1）Landsat TM影像包含7个波段，波段1-5（可见光部分）和波段7（短波红外）的空间分辨率为30米，波段6（热红外波段）的空间分辨率为120米。2）ASTGTM_N35E

7、103条带号：103 行编号：35 中心经度：103.5 中心纬度：35.5 分辨率30米数据投影为UTM/WGS844. 研究区域：ASTGTM_N35E103研究区域是甘肃省兰州市榆中县地区，中心经纬度是35.5N，103.5E。榆中县隶属于甘肃省中部，榆中县地势南高北低，中部凹，呈马鞍形。南部为石质高寒山区，马寒山最高峰海拔3670.3米。北部为黄土丘陵区，最高峰吕家岘海拔2495米。南北两山之间是川塬丘陵沟壑区，海拔1500-2000米，地形由西南、东南和东北三面向西北倾斜。榆中县属于温带半干旱气候，每年均气温6.7，降水量400毫米，无霜期120天。榆中县有耕地105.74万亩，占全

8、市总耕地面积的33.4%, 其中有效灌溉面积29.14万亩。榆中县内已探明的矿藏15种，有中小型矿床，矿点和矿化点48处。矿产种类繁多，已探明的有48种。操作步骤1. 从下载的ASTER DEM中选取512512子区数据，显示分析研究区地形地貌特征a. 打开地理空间数据云网站，注册用户，在高级检索里面，搜索数据集（DEM数字高程数据GDEMDEM 30M分辨率数字高程数据），搜索甘肃省兰州市榆中县影像，ASTGTM_N35E103b. 打开ENVI软件，利用Resize Data（Spatial/Spectral）工具从多波段影像中选取512512子区并保存。打开Resize Data Inp

9、ut File点击Spatial Subset出现Select Spatial Subset对话窗口，继续点击Image,出现带图的对话窗口，将研究范围大小改为512*512（边长应控制在5121024之间）拖动正方形至有山地起伏的地方，选择保存路径保存为512。2. 完成高通滤波（High Pass），拉普拉斯变换（Laplacian），定向增强（Directional），高斯高通滤波（Gaussian High Pass），Sobel边缘增强（Sobel），Roberts算子增强（Roberts）增强运算，比较变换前后影像中地形/地貌信息上的差异。a. 点击工具栏Filter-Convol

10、utions and Morphology（卷积与形态），打开Convolutions and Morphology Tool窗口，Convolutions可以选择10中卷积方法，分别为High Pass（高通滤波）， Low Pass（低通滤波器）， Laplacian（拉普拉斯变换）， Directional（定向增强）， Gaussian High Pass（高斯高通滤波） ，Gaussian Low Pass（高斯低通滤波器）， Median（中值滤波变换） ，Sobel（Sobel边缘增强）， Roberts（Roberts算子增强）， User Defined（用户自定义）b. K

11、ernel Size可以改变卷积核大小，（以奇数来表示如3*3,5*5，有的卷积是不能改变大小的，如Sobel, Roberts，如需要使用非正方形，应选择Options-Square Kernel）c. Image Add Back可以改变加回值，将原始图像中的一部分“加回”到卷积滤波结果图像上，有助于保持图像的连续性，该方法经常用于图像锐化。“加回值”是原始图像在结果输出图像中所占百分比。d. Editable Kernel卷积核中各项的值。在文本框中双击鼠标可以进行编辑，选择File-Save Kernel或Restore Kernel，可以把卷积核保存为文件（.ker）或者打开一个卷积

12、核文件e. 这里我们分别做要求的前4种卷积运算的3中形式，分别为Kernel Size为5*5，加回值0；Kernel Size为7*7，加回值0；Kernel Size为7*7，加回值30。和Sobel, Roberts增强的加回值分别为0,10,30的运算形式。其他默认。f. 保存相关文件，并进行比对。5. 纹理运算，比较变换前后影像中地形/地貌信息上的差异：l 一阶纹理运算l 二阶纹理运算一阶纹理运算：a. 点击工具栏File-Occurrence Measures，打开Occurrence Texture Parameters窗口b. Textures to Compute(结构计算)

13、里有5个不同的基于概率统计的纹理滤波：Data Range数据范围，Mean平均值，Variance方差，Entropy信息熵和Skewness偏斜c. Processing Window为处理窗口的大小，此处默认d. Choose选择要保存的途径与文件名Occurrencee. 设立6个Views，在每个窗口分别打开原图与5个纹理滤波处理后的结果二阶纹理运算：a. 点击工具栏File-Co-occurrence Measures，打开Co-Occurrence Texture Parameters窗口b. Textures to Compute(结构计算)里有8个不同的基于二阶矩阵的纹理滤波

14、：均值Mean，方差Variance，协同性Homogeneity，对比度Contrast，相异性Dissimilarity，信息熵Entropy，二阶矩Second Moment和相关性Correlationc. 在“Rows”(Y)和“Cols”(X)文本框中，键入处理窗口的大小d. 键入X,Y变换值（Co-occurrence Shift），用于计算机二阶概率矩阵e. Greyscale quantization levels选择灰度量化级别：None,64,32或者16f. 选择输出路径及文件名Co-occurrence。g. 设置6个Views，分别加载除却与一阶纹理运算重复的剩余的

15、5个特征。3. 周期噪声去除：利用傅立叶变换对影像tm_1.img所有7个波段进行处理，选择恰当的参数消除影像上的噪声信号a. 打开ENVI Classic，点击工具栏File-Open Image File打开TM_1影像b. 在Available Bands List中加载TM_1影像的Layer_1，右击Layer_1，点击Load Band to New Display，在新窗口中加载Layer_1c. 在总窗口工具栏Filer-FFT Filtering-Forward FFT，选择TM_1.img，点击OK。选择输出文件路径与名称fft，其他默认。在Available Bands

16、List中加载fft影像的Layer_1d. 在总窗口工具栏Filer-FFT Filtering-Filter Definition（筛选定义），在随后的窗口中选择第一个，点击OK。在Filter Definition窗口中，在Filter_ Type下选择User Defined Cut用户自定义，其他默认e. 在fft Layer_1的子窗口工具栏中，点击Overlay-Annotation，打开Annotation Text窗口，在Object下选择Polygon，打开镜像功能，在Color中定义滤波器的颜色，随意定义只要有颜色凸显。在fft Layer中绘制滤波器。滤波器要求尽量没有

17、空隙，上下又不交叉。（如有空隙则消除噪声不干净，还会有类似条带状的纹路。如有交叉重叠，则像元会有重叠，消除噪声后的结果明显混乱）。右键可以镜面对称，右击绿色圈圈可以上下对称。做好滤波器后，在子窗口工具栏File-Save Annotation保存文件为b1_fft好的滤波器绘制如图下f. 回到Filter Definition窗口，在Ann File下导入刚才保存的文件，在Choose下保存定义好的滤波器名为B1。Apply应用。此时B1滤波器已出现在Available Bands List中点击可以加载查看。g. 在总窗口工具栏Filer-FFT Filtering-Inverse FFT，

18、在Inverse FFT Input File中选择变换后的影像fft，在Spectral Subset中选择与之对应的波段Layer_1，点击OK。在随后出现的窗口选择滤波器B1，点击OK。选择输出径Band_1。此时波段1处理后的结果Band_1已出现在Available Bands List中点击可以加载查看。h. 重复以上步骤，分别处理其他波段，得到Band_1,2,3,4,5,6,7。i. 在总窗口工具栏第二个选项中打开Layer stacking，打开窗口，点击Input File，选中Band1,2,3,4,5,6,7。在Boarder File中可以给波段排序。在Choose中

19、选择输出路径并保存文件为jieguo。到此为止便生成了ENVI的标准格式。j. 在Available Bands List中点击RGB Color真彩色。以432波段分别加载原始影像与处理后的影像。结果与分析一 截得ASTGTM_N35E103的512子区域由图可见山地高耸区域呈深黑，平坦地区呈高亮，山脊为树枝状二 卷积增强运算部分1. 高通滤波（High Pass）处理后的结果a. 用高通滤波处理后，整体纹理粗糙，明暗对比度变小，呈现不明显的浮雕效果，坡地凸出。容易区分沟谷区域。b. 比较high7_0与high5_0即不加回值的两个。可知窗口变大时，边缘锐化程度明显,浮雕效果与立体效果更加

20、明显。c. 比较high7_0与high7_30即只有加回值不同的两个。可知加回值变大，图像边界的明暗对比度增强（主要由图像右上角得知）d. 总结：高通滤波使影像呈现浮雕般立体感，立体感的强弱取决于卷积窗口的大小，窗口越大，立体感越明显。高通滤波在保持图像高频信息的同时，消除了图像中的低频成分。它可以用来增强纹理、边缘等信息。输出的影像和原影像相比，增强了纹理和边缘信息，有明暗交错的效果。随着卷积核大小的增加，图像的边缘被凸显的越来越明显。特征被进一步放大。2. 拉普拉斯变换（Laplacian）处理后的结果a. 拉普拉斯检测的是变化的频率，是一个二阶的概念。白色斑块处代表地形起伏剧烈的地方。

21、主要集中在山地地区。b. 不加回值改变窗口大小，窗口越大，浮雕效果越明显c. 加回值窗口大小相同，图像边界明暗对比度增强，对沟谷表现能力增强d. 总结：拉普拉斯滤波调大影像加回值后，在影像中能更清晰的反映地势较低的区域。影像整体对比度增强，在影像中能更清晰的体现出地势变化。拉普拉斯滤波是边缘增强滤波，它的运行不用考虑边缘的方向。输出的影像和高通滤波处理过的影像相似，但是更加平滑，增强了纹理和边缘信息，河流、山脊、沟谷等边缘特征被凸显。3. 定向增强（Directional）处理后的结果a. 影像立体感较强，呈现较强的浮雕效果。山脊尤为凸出，能够清晰地反映出山体边缘，沟谷及褶皱。b. 同样遵循加

22、回值，窗口改变的带来的规律。对比加回值，窗口都一致，方向分别为0度，45度，90度，处理后的结果。c. 可以发现随度数的增加，浮雕效果减弱，明暗对比减小。d. 0度时，更能凸出山脊，山谷及褶皱。东边高亮，西边暗一些。（可以由Editable Kernel中的+1-1得知，+1高亮-1变暗）e. 左下角90度时，南边高亮，北边会暗一些f. 总结：方向滤波器可根据不同角度的定向增强运算使影像不同方向体现不同的明暗效果，使影像立体感更强。在不同的明暗效果中可分析河流、沟谷、山体褶皱的具体形态。方向滤波器是边缘增强滤波，它有选择性地增强有特定方向成分的图像特征。输出影像呈现浮雕效果，立体感比较强。4.

23、 高斯高通滤波（Gaussian High Pass）处理后的结果a. 高斯高通滤波器通过一个指定大小的高斯卷积函数对图像进行滤波。高斯高通滤波处理对于突出体现影像中高频区域有较好的效果。影像整体偏亮，影像纹理粗糙，河流、沟谷等均不明显b. 加回值不变，窗口变大，浮雕效果变得明显。c. 在7*7窗口上加30%回值，影像边缘更平滑，地势起伏变化更加流畅，与原图较为相似。5. Sobel边缘增强（Sobel）处理后的结果a. 滤波器是非线性边缘增强滤波器，它是使用 Sobel 函数的近似值的特例，滤波器的大小不能更改，也无法对卷积核进行编辑。b. Sobel加回值为0时, 此操作后影像与原DEM影

24、像差异较大，影像纹理粗糙，呈现不规则的线状。影像突出了地物的边界信息，如河流沟谷、山体边缘等，但对于影像地势体现不明显，影像立体感较弱。c. Sobel加回值为30时，此操作后影像中地物边界信息表现更为明显，颜色对比度更加明显，且边缘锐化增强。对于山体起伏表现较弱。6. Roberts算子增强（Roberts）处理后的结果与Sobel对比a. Roberts加回值为0时，地物边缘同样得到锐化。但影像中地势明暗对比度较弱，因此对于山地区域地势变化难以分辨。且操作后地物边缘清晰度较差。b. Roberts加回值为30时，与影像加回值0%时相比，此操作后影像立体感增强，地物边缘凸显且更加清晰。在影像

25、中能较好的分辨沟谷。且影像对比度增强，能较好的体现地势变化。c. Roberts加回值为0与Sobel加回值为0很像，但前者边界对比度稍微强一些。d. robert滤波器：地物边界得到了锐化，随着影像加回值的增加，输出影像的平滑度增加，影像越来越模糊，细节信息减少。二 纹理运算部分1.一阶纹理运算a. Data Range数据范围地势起伏的地方灰度变化大，所以高亮显示；平坦的地方，灰度变化小，所以区域显示黑暗。立体效果不明显，边界明显。影像中河流沟谷凸显，山体边缘清晰，山体褶皱呈现线状，影像整体立体感较差。b. Mean均值与原图类似，但是Mean放大后可以发现，它在边界处相对比较平滑，细节部

26、分消失c. Variance方差与Data Range类似，平地区域地势平坦，值方差较小，故在影像中呈现暗黑色；沟谷地区坡度变化较大，在影像中高亮显示。但可以看出它整体相对来说暗一些，且边缘更为突出d. Entropy信息熵整体偏亮，山地褶皱地区DEM值变化信息复杂，则此区域熵值大，在图像中高亮显示；而在沟谷、平原区域，地势平坦，地物信息变化小，则此区域显示暗黑。e. Skewness偏斜整体色调为灰色，偏度表示地物的倾斜程度。值大于零表示地物向左偏斜，值小于零表示地物向右偏斜，值等于零表示地物对称分布。图中高亮部分代表地形倾斜，不对称。2. 二阶纹理运算a.协同性Homogeneity:表示影像特定区域范围内灰度值的相似性。由于在地势平坦，则区域同质性大，故高亮显示；在沟谷山地区域内地势崎岖，同质性小，故显示暗黑。b.对比度Contrast: 体现图像纹理沟深，纹理相对粗糙