【ENVI入门系列】21.波段运算与波谱运算工具

1
分享 2016-06-14
版权声明:本教程涉及到的数据仅供练习使用,禁止用于商业用途。
目录
波段运算与波谱运算工具
1. 概述
2. 波段运算用法示例
3. 波段运算满足条件
4. 波段运算的IDL知识
4.1 注意数据类型
4.2 数据类型的动态变换
4.3 充分利用IDL的数组运算符
4.4 运算符操作顺序
4.5 调用IDL函数时注意事项
5. 波段运算典型公式
6. 调用IDL用户函数
6.1 编写函数
6.2 编译函数
6.3 使用函数
7. 波谱运算(Spectral Math)


1. 概述
ENVI Band Math是一个灵活的图像处理工具,其中许多功能是无法在任何其它的图像处理软件中获得的。由于每个用户都有独特的需求,利用此工具,用户自己可以定义处理算法,应用到ENVI打开的波段或整个图像中,用户可以根据需要自定义简单或复杂的处理程序。例如,可以对图像进行简单加、减、乘、除运算,或使用IDL编写更复杂的处理运算功能。
波段运算实质是对每个像素点对应的像素值进行数学运算。如下图为一个简单波段运算的示意图,运算表达式是三个变量相加,每一个变量对应于一个图像数据,对这三个图像数据求和并输出结果图像。表达式中的每个变量不仅可以对应于单一波段,也可以是一个多波段的栅格文件。例如:在表达式b1+b2+b3中,如果b1是一个多波段图像文件,b2、b3为单一波段,则结果为b1所对应的文件的所有波段分别和b2、b3进行求和。

图:波段运算示意图

2. 波段运算满足条件
使用波段运算需要满足5个基本条件:
一、必须符合IDL语言书写波段运算表达式
所定义的处理算法或波段运算表达式必须满足IDL语法。不过,书写简单的波段运算表达式无须具备IDL的基本知识,但是如果所感兴趣的处理需要书写复杂的表达式,建议学习用于波段运算的IDL知识。
二、所有输入波段必须具有相同的空间大小
由于波段运算表达式是根据pixel-for-pixel原理作用于波段的,因此输入波段在行列数和像元大小必须相同。对于有地理坐标的数据,如果覆盖区域一样,但是由于像元大小不一样使得行列数不一致,在进行波段运算前,可以使用/Raster Management/Layer Stacking功能对图像进行调整。
三、表达式中的所有变量都必须用Bn(或bn)命名
表达式中代表输入波段的变量必须以字母"b"或"B"开头,后跟5位以内的数字。例如:对3个波段进行求和运算的有效表达式可以用以下3种方式书写:
  • b1 + b2 + b3
  • B1 + B11 + B111
  • B1 + b2 + B333


四、结果波段必须与输入波段的空间大小相同
波段运算表达式所生成的结果必须在行列数方面与输入波段相同。例如,如果输入表达式为MAX(b1),将不能生成正确结果,因为表达式输出值为一个数,与输入波段的行列数不一致。
注:MAX函数功能为求数组最大值。
五、调用IDL编写的自定义函数时
波段运算工具可以调用IDL编写的Function,当函数为源码文件(.pro)时,必须启动ENVI+IDL才能调用;如果函数编译为了sav文件,可以将sav文件放到如下路径,重启ENVI即可调用。
  • ENVI 4.x: C:\Program Files\ITT\IDL\IDL80\products\envi48\save_add
  • ENVI Classic: C:\Program Files\Exelis\ENVI51\classic\save_add
  • ENVI 5.x: C:\Program Files\Exelis\ENVI51\extensions



3. 波段运算用法示例
下面以求一个图像数据三个波段的和为例,学习Band Math工具的使用。示例数据为ENVI自带数据,位于"…\21.波段运算与波谱运算工具\数据\can_tmr.img"。
(1) 启动ENVI,选择菜单File > Open,打开数据"can_tmr.img";
(2) 启动Band Math工具,路径为Toolbox/Band Ratio/Band Math;
(3) 在Band Math面板,在Enter an expression文本框中输入运算表达式:b1+b2+b3,点击Add to List按钮,将表达式添加到Previous Band Math Expression列表中;
注:如果表达式存在语法错误,将不能被添加到列表中,如下图所示。

图:Band Math面板
(4) 在Band Math面板中,选中添加的"b1+b2+b3",点击OK按钮,打开Variables to Bands Pairings对话框(如图),为运算表达式中各个变量赋图像文件或者图像波段;

图:Variables to Bands Pairings对话框
(5) 在Variables to Bands Pairings对话框中,Variables used in expression列表框中选择变量b1,单击Available Bands List中的"TM Band 1 (0.4850)"。然后用同样的方法为b2和b3指定为"TM Band 2"和"TM Band 3";
注1:本例中指定变量为波段,同样可以指定为文件。方法如下:选中变量后,点击按钮Map Variable to Input file,可以为变量指定一个多波段图像文件
注2:一旦第一个波段或文件被选中,只有那些具有相同行列数的波段被显示在波段列表中。
(6) 单击Choose按钮,选择文件名及路径保存结果,单击OK按钮执行运算。
(7) 此时可以将输入和输出文件加载到视图中,然后点击工具栏
图标获取当前鼠标位置的像元值,如图所示。

图:鼠标取值Cursor Value面板
注:本例中输入文件的数据类型为Byte型,范围0~255,所以使用b1+b2+b3有可能会越界。所以,更标准的公式应该为fix(b1)+b2+b3。本例只是简单演示用法。后续章节会详细介绍波段运算的注意事项,如数据类型的转换等。
在Band Math对话框中,以下是其他按钮的功能说明:
  • 单击Save按钮可以将列表中的运算表达式保存为外部文件(.exp);
  • 单击Restore按钮可以将外部运算表达式文件导入;
  • Clear按钮可以清除列表中的所有运算表达式;
  • Delete按钮可以删除选择的运算表达式。



4. 波段运算的IDL知识
波段运算的强大功能是由IDL的功能、速度和灵活性所提供的。但是要熟练使用波段运算功能,并不需要成为一个熟悉IDL编程的专家。下面的知识可以帮助熟练使用波段运算功能并避免一些经常出现的问题。

4.1 注意数据类型
IDL中的数学运算与简单的使用计算器进行运算是有一定差别的。要重视输入波段的数据类型和表达式中所应用的常数。每种数据类型——尤其是非浮点型的整型数据都包含一个有限的数据范围。例如:8-bit字节型数据表示的值仅为0-255,如果对波段求和(b1 + b2)并且其值大于255,那么得到的结果将与期望值不等。当结果值大于某个数据类型所能容纳的值的范围时,该值将会溢出(overflow)并从头开始计算。例如,将8-bit字节型数据250和10求和,结果为4。
IDL> a = 250B
IDL> b = 10B
IDL> help, a, b
A BYTE = 250
B BYTE = 10
IDL> print, a + b
4
类似的情况经常会在波段运算中遇到,因为遥感图像通常会被存储为8-bit字节型或16-bit整型。要避免数据溢出,可以使用IDL中的一种数据类型转换功能(参见下表)对输入波段的数据类型进行转换。例如:在对8-bit字节型整型图像波段求和时(结果有大于255),如果使用IDL函数 FIX() 将数据类型转换为整型,就可以得到正确的结果。
fix(b1)+ b2
你可能会有这样的想法:既然浮点型或双精度数据可以表示所有的数据值,为什么不在所有的计算中都使用浮点型数据呢?这是因为一个数据所能表现的动态数据范围越大,它占用的磁盘空间越多。例如:字节型数据的一个像元仅占用1个字节;整型数据的一个像元占用2个字节;浮点型数据的一个像元占用4个字节。浮点型结果将比整型结果多占用一倍的磁盘空间。关于IDL数据类型的占用磁盘空间和数据范围的详细介绍,参考下表。
表:数据类型及说明


4.2 数据类型的动态变换
一些数字可以使用几种不同的数据类型表达出来,IDL制定了一些默认规则对这些数据进行解译。因此IDL的数据类型是可以进行动态变换的,也就是说IDL能够将表达式中的数据类型提升为它在表达式中所遇到的最高数据类型。例如,一个整型数字,即使它在8-bit字节型的动态范围,也常被解译为16-bit整型数据。如果想为一幅8-bit字节型数据图像加5,并且使用如下的波段运算表达式:
B1 + 5
数据5将被解译为16-bit整型数据,因此波段运算结果将被提升为16-bit整型数据图像(占用8-bit字节型图像的两倍磁盘空间)。如果想保持结果为字节型图像,可以使用数据类型计算函数byte( ):
b1 + byte(5)
或使用IDL中将16-bit整型数据转换为8-bit字节型数据的缩写:
b1 + 5B
在数据后紧跟一个字母B表示将该数据解译为字节型数据。如果在波段运算表达式中经常使用常数,这些类似的缩写是很有用的。所有数据类型的缩写见上表。

4.3 充分利用IDL的数组运算符
IDL的数组运算符使用方便且功能强大。它们可以对图像中的每一个像元进行单独检验和处理,因而避免了FOR循环的使用(不允许在波段运算中使用)。数组运算符包含关系运算符(LT、LE、EQ、NE、GE、GT)、Boolean运算符(AND、OR、NOT、XOR)和最小值、最大值运算符(<、>)。这些特殊的运算符对图像中的每个像元同时进行处理,并将结果返还到与输入图像具有相同维数的图像中。
例如:要找出所有负值像元并用值-999代替它们,可以使用如下的波段运算表达式:
(b1 lt 0)*(-999)+(b1 ge 0)* b1
关系运算符对真值(关系成立)返回值为1,对假值(关系不成立)返回值为0。系统读取表达式(b1 lt 0)部分后将返还一个与b1维数相同的数组,其中b1值为负的区域返回值为1;其他部分返回值为0,因此在乘以替换值-999时,相当于只对那些满足条件的像元有影响。第二个关系运算符(b1 ge 0)是对第一个的补充——找出那些值非负的像元,乘以它们的初始值,然后再加入替换值后的数组中。类似的使用数组运算符的表达式为波段运算提供了很强的灵活性。下表中描述了IDL数组操作函数。
表:IDL数组操作函数和运算符



4.4 运算符操作顺序
在波段运算过程中,是根据数学运算符的优先级对表达式进行处理,而不是根据运算符的出现顺序。使用圆括号可以更改操作顺序,系统最先对嵌套在表达式最内层的部分进行操作。IDL运算符的优先级顺序列在下表中。具有相同优先级的运算符根据它们在表达式中出现的顺序进行操作。例如,考虑如下表达式(用常数代替波段):
  • 5 + 3 * 2 —— 求得的值为11,因为乘号运算符的优先级高
  • (5 + 3)* 2 —— 求得的值为16,因为圆括号改变了操作顺序


将优先级的顺序与数据类型的动态变换结合起来时,如果操作不当,也将改变表达式的运算结果。要确保将表达式中的数据提升为适当的数据类型,从而避免数据的溢出或在处理整型除法时出现错误。例如,考虑如下示例:
float(5)+ 10 / 3
所有的常数都为整型,但 float() 函数将结果提升为浮点型数据,由于除号的优先级高于加号,因此先以整型数据进行除法运算,将结果与被提升为浮点型数据的5相加得到一个浮点型结果8.0,而不是所期望的结果8.3。
5 + 10 / float(3)
如果将数据类型转换函数移到除法运算中,将得到期望的结果8.3。
表:运算符优先级


4.5 调用IDL函数时注意事项
如同其他ENVI程序一样,波段运算处理也是分块进行的。如果被处理的图像大于在参数设置中被指定的碎片(tile)尺寸,图像将被分解为更小的部分,系统对每一部分进行单独处理,然后再重新组合起来。当使用的IDL函数同时需要调用所有图像数据时,由于波段运算表达式是对每一部分数据进行单独处理的,这种处理方法将会产生问题。例如,在使用求取数组中的最大值的IDL函数MAX()时:
b1 / max(b1)
如果波段运算是分块进行的,则每一个部分除以的值是该部分的最大值,而不是整个波段的最大值。如果运行这个运算式发现波段运算结果中有较宽的水平条带,那很有可能是由于分块处理造成的,因为图像是水平分块的。
其他需要注意的IDL函数还包括:MAX、MIN、MEAN、MEDIAN、STDDEV、VARIANCE和TOTAL等。在多数情况下,使用BYTSCL函数也比较困难。

5. 波段运算典型公式
一、避免整型数据除法
当对整型数据波段进行除法运算时,运算结果不是被向上或向下取整,而是直接被简单地舍去(小数点后面的数据被舍弃)。要避免这种情况发生,通常将数据类型转换为浮点型。
b1 / float(b2)
如果想将除法数据结果保持为整型,最好先将数据转换为浮点型进行除法运算,然后再将结果转换为所需的数据类型。例如:如果输入波段为8-bit字节型,想将结果取整并存储为16-bit整型数据,使用下面的表达式:
fix(ceil(b1 / float(b2)))
二、避免整型运算溢出
整型数据包含一个动态的数据范围。如果波段运算将生成的数据相当大或相当小,无法以输入波段的数据类型表示出来,要注意提升相应的数据类型。例如:如果示例表达式中的波段b1和b2为8-bit字节型数据,生成结果的最大值可能为(256*256)=65,025。由于字节型数据所能表示的最大值为255,因此结果的数据类型只有被提升为16-bit无符号整型才能返回正确的值,否则,大于255的值将溢出,并记录一个错误的值。因此使用下面表达式避免数据范围溢出:
uint(b1) * b2
三、生成混合图像
波段运算为多幅图像的混合提供了简单的方法。例如:如果b1和b2为8-bit字节型数据,下面的表达式将生成一幅新的8-bit字节图像,b2所占权重为0.8,b1所占权重为0.2。
byte(round((0.2 * b1) + (0.8 * b2)))
四、使用数组运算符对图像进行选择性更改
波段运算为图像的选择性更改和来自多幅图像的数据结合提供了简单的方法。在下面的示例中,把两幅图像结合起来进行处理,从而从图像中消除云的影响。在图像b1中,像元值大于200的像元被认为是云,希望用图像b2中的相应像元对它们进行替换。
(b1 gt 200)* b2 + (b1 le 200)* b1
用类似的运算表达式,可以将一幅图像的黑色背景变成白色背景:
(b1 eq 0)*255 + (b1 gt 0)*b1
下面的示例是一个较为复杂的表达式。该表达式使用几个标准来生成一幅二进制掩膜图像,用于识别主要为云的像元。该算法可以应用于经过定标的AVHRR日间图像中生成云的掩膜图像。在该表达式中,b4(热红外波段)值必须为负,或b2(反射波段)值必须大于0.65并且b3和b4(中红外和热红外波段)的差值必须大于15度。由于关系运算符为真值(关系成立)返回1值,因此生成的掩膜图像在有云处值为1,在其他区域值为0。
(b4 lt 0) or (b2 gt 0.65) and (b3–b4) gt 15
五、最小值和最大值运算符的使用
最小值和最大值运算符也是数组的基础运算符,但与关系运算符或Boolean运算符不同的是:它们不返还真值或假值,而返还实际的最小值和最大值。在下面的例子中,对于图像中的每一个像元,0、b2或b3中的最大值将被加到b1中,该表达式确保加到b1中的值始终为正。
b1 +(0 > b2 > b3)
在下面的例子中,最小值和最大值运算符的同时运用使b1中的值被限制在0和1之间,最后得到的结果在[0,1]范围内。
0 > b1 < 1
有时候需要计算几年内的数据平均值(如NDVI),如果某数据的值为0则不参加计算,如果3个通道都为0,则赋值为0,比如某点b1=4;b2=6;b3=0;那么平均值ave = (b1+b2+b3)/(1+1),则可用以下运算表达式:
(b1>0+b2>0+b3>0) / (((b1 ge 0) + (b2 ge 0)+(b3 ge 0)) >1)
六、利用波段运算修改NaN
NaN为Not a Number的缩写,在遥感图像中属于异常值。很多用户有修改NaN的需求,比如把0值修改为NaN,或把NaN修改为0值等。由于波段运算公式较为复杂,现归纳如下。
  • 修改0值为NaN

       1)float(b1)*b1/b1
  • 修改特定值(250)为NaN

       2)b1*float(b1 ne 250)/(b1 ne 250)
  • 修改NaN为特定值(-999):

       3)finite(b1, /nan)*(-999) or (~finite(b1, /nan))*b1
  • 修改NaN为0值(先按上面方法修改为-999或其他图像中不存在的值)

       4)(b1 ne -999)*b1


6. 调用IDL用户函数
ENVI 提供对 IDL 程序的访问的功能,可以使用内置的IDL 函数或者用户自定义IDL函数。这些函数要求它们接受一个或多个图像阵列作为输入,并且输出一个与输入波段具有相同行列的单波段二维数组作为计算结果。如下为一个自定义函数的基本格式:

FUNCTION bm_func, b1, [b2,..., bn, parameters and keywords]
processing steps
RETURN, result
END



下面以一个简单的例子介绍用IDL自定义函数,后在Band Math中使用这个函数,自定义函数实现的功能是计算一个比值(b1+b2)/(b1-b2),并且检查分母为0的情况。

6.1 编写函数
用记事本或IDL工作台新建文件,编写以下代码:

;定义两个变量和一个关键字
FUNCTION bm_ratio, b1, b2, check=check
;计算分母
den = FLOAT(b1) - b2
IF (KEYWORD_SET(check)) THEN $
ptr = WHERE(den EQ 0, count) $
ELSE $
;如果设置了check关键字,检查分母为0情况
count = 0
;如果分母为0,临时则将分母赋值1.0
IF (count GT 0) THEN den[ptr] = 1.0
;计算比值
result = (FLOAT(b1) + b2) / den
;分母为0时,直接将结果返回0.0
IF (count GT 0) THEN result[ptr] = 0.0
RETURN, result
END



将源码文件保存为"bm_ratio.pro"。

6.2 编译函数
有三种方式编译这个自定义函数:
(1) 将bm_ratio.pro文件拷贝到安装路径的Extensions或save_add 目录下,启动ENVI + IDL模式,自动将bm_ratio.pro编译。
(2) 启动 ENVI + IDL模式,在IDL编辑器中打开bm_ratio.pro源码文件,点击IDL工作台工具栏的编译按钮,然后就可以在ENVI中使用bm_ratio函数了。
(3) 在IDL中,通过save命令,将pro源码文件编译为sav文件,然后拷贝到ENVI安装路径的Extensions或save_add文件夹下,启动ENVI即可使用。

6.3 使用函数
(1) 打开一个多光谱文件;
(2) 与第2章节方法一样,打开波段运算工具,输入如下表达式:
bm_ratio(b1, b2) ; 不执行检查分母为0的情况

bm_ratio(b1, b2,/check) ; 执行检查分母为0的情况

图:使用IDL用户函数
(3) 其他操作过程与第2章节一样,为变量选定波段或文件后,执行即可。

7. 波谱运算(Spectral Math)
ENVI Spectral Math是一种灵活的波谱处理工具,可以用数学表达式或IDL程序对波谱曲线(以及选择的多波段图像)进行处理。波谱曲线可以来自一幅多波段图像的Z剖面、波谱库或ASCII文件。
如图所示为波谱运算的简单示意图——求三个波谱曲线的和。在表达式s1+s2+s3中(波谱运算中的变量是以s开头),可以分别给s1、s2、s3指定为一条波谱曲线,得到的结果是一条波谱曲线(x值与s1、s2、s3一样,y值是三者之和);也可以s1是一个多波段图像文件(其实是每个像素点的Z-剖面),s2和s3分别是两条波谱曲线,得到的结果是一个与输入的多波段图像一样波段数和行列数的图像。

图:波谱运算示意图
下面以求表达式 (s1+s2+s3)/3 为例介绍Spectral Math工具的使用,输入数据源为ENVI自带的波谱库文件。
(1) 选择菜单Display > Spectral Library Viewer,启动波谱库浏览器;
(2) 在Spectral Library Viewer面板左侧,自动加载了ENVI自动波谱库文件。这里选择veg_lib文件夹里的veg_1dry.sli波谱库文件;
(3) 左键在波谱列表中单击,光谱曲线将显示在面板右侧的视图中,任意选择三个光谱曲线,结果如下图所示,可以点击右侧中间的三角箭头浏览已选的波谱列表;
(4) 启动Spectral Math工具,路径为Toolbox/Spectral/Spectral Math;
(5) 在弹出的Spectral Math面板中,输入表达式(s1+s2+s3)/3,点击Add to List将公式添加到上面的列表中,选中列表中的公式,点击OK;
注:由于波谱库中的光谱数据类型为浮点型,所以不需要进行数据类型转换。如果输入光谱数据类型为字节型、整型等,需要进行数据类型转换。

图:Spectral Library Viewer面板
(6) 在弹出的Variables to Spectra Pairings面板中,为s1、s2、s3指定波谱曲线。Output Result to可以选择Same Window或New Window,这里选择New Window,将结果输出到新窗口中。可通过鼠标左键将输入三条曲线拖拽到新的窗口中,结果如图所示,红绿蓝三条曲线为输入波谱,紫色为输出波谱。

图:Variables to Spectra Pairings面板

图:波谱运算结果
注:波谱运算工具与波段运算一样,可以调用IDL编写的用户函数。编写和调用方法与波段运算工具一致,这里不再赘述。

教材下载:http://pan.baidu.com/s/1jGgM1cq
数据下载:http://pan.baidu.com/s/1sjJNQvb
视频下载:http://pan.baidu.com/s/1eQu1KkU

文章来源:http://blog.sina.com.cn/s/blog_764b1e9d0102v58s.html

0 个评论

要回复文章请先登录注册