基于IDL的遥感影像波段运算

1.IDL语言特征
在IDL中，矩阵是按照列的方式处理的，即矩阵是以n列、m行的方式表示的，这一点与Fortran语言一样，但与C语言（按m行、n列的维度来标记）表示矩阵的方式不同。只有熟悉IDL的数据存储、处理的特点，才能方便地进行分析。以300列、200行的遥感数据band1为例，其数据的记录特征如下：

像元坐标
亮度矩阵

1,1 2,1 …… 300,1 92 94 …… 90
1,2 2,2 …… 300,2 96 92 …… 91
…… …… …… …… ……
1,200 2,200 …… 300,200 98 101 …… 90






IDL图形坐标 band1 DN值坐标
图1.IDL中像元坐标与亮度矩阵的对应关系
（1）数据按列存储；最先存储的是第一列的数据（图象第一列的像元亮度值），接着是第二列，直至最后一列。
（show:print,data[0,*]）
（2）采用顺序显示数据的方式。将第一行的亮度值按顺序从左到右以行的方式显示或存储，不管band1矩阵有多少行、多少列，总是先显示、存储完第一行，然后是下一行。由于IDL的默认显示宽度是80列，所以在输出窗口中要用几行才能显示完矩阵的第一行亮度值。然后另起一行显示亮度矩阵的下一行值，直到最后一行。这一点对于图象显示和分析相当重要。
复习：IDL矩阵的显示方式，data
（show:print,data[*,0]）
（3）默认情况下，图象数组的第一行和第一个元素显示在屏幕的左下角。后面的行将从下往上显示。在图形窗口显示图象时，图象的每一个元素在屏幕中显示为一个像素[1]。在IDL中图象是从下往上画的。即把band1亮度矩阵的第一行值显示在图象最下端的一行像元中，把band1亮度矩阵的最后一行值显示在图象最上端的一行像元中。这是图象显示最重要的特征。
区别两幅图象的显示效果：（tvscl,data）与（tvscl,data,300，0，order=1）
（4）矩阵转置与图象的关系
将矩阵band1转置，
(print,(transpose(data))[*,0])，打印转置矩阵的第一行，它在band1中是哪行或哪列？
(print,transpose(data)
tvscl,data,200,0,order=1
结论：按行显示的数据DN，按列显示；按列显示的DN数据，按行显示
2.遥感数据的处理与分析
31通用的栅格数据存储格式
BSQ ：波段顺序存储格式。每行数据后面紧接着同一波谱波段的下一行数据。这种格式最适于对单个波谱波段中任何部分的空间（X，Y）存取。每个tile是单个波段的一个空间子集。
BIP ：波段按像元交叉格式。图像按顺序存储第1个像元所有的波段，接着是第2个像元的所有波段，然后是第3个像元的所有波段，等等，交叉存取直到像元总数为止。这种格式为图像数据波谱（Z） 的存取提供最佳性能。每个tile是一幅图象中所有波段行的所有像素。
BIL：波段按行交叉格式。按 BIL 格式存储的图像先存储第一个波段的第一行，接着是第二个波段的第一行，然后是第三个波段的第一行，交叉存取直到波段总数为止。每个波段随后的行按照类似的方式交叉存取。这种格式提供了空间和波谱处理之间一种折衷方式，它是大多数 ENVI 处理任务中所推荐的文件格式。每个tile是一幅图像所有波段的一行。
32栅格数据输入
以BSQ格式的图象为例。假设有一幅BSQ格式的7波段TM影象“2003_sw.img”，需要计算各波段亮度值的基本统计信息（包括：最大值、最小值、均值、均方差、倾斜度和平坦度）。
遥感数据其实就是像元的亮度值（也叫做DN值），这与它有多少个波段、是多光谱数据还是高光谱数据都没有关系。这里分析遥感数据的目的，是为了帮助理解、纠正易于混淆的栅格数据存储方法和数据分析方法。影象“subset_2003.img”的大小为300×200（即：lines:200，Samples：300）。以band1波段为例，在IDL的图形显示系统中，左上角坐标为（1，1）。
利用程序把band1的亮度值读入IDL环境中，构造成一个亮度矩阵（见图1）。在IDL的图形系统中，图象左上角的坐标为（1，1），右下角（X,Y）坐标值最大。在图1中，图象band1的像元坐标左上角为（1，1），右下角为（300，200），坐标对中的X为列号，Y为行号，像元数量为300列200行。与像元一一对应，亮度矩阵是band1的亮度值，每个像元一个值，左上角为92，右下角为90。亮度矩阵中的第一行值，与图象中的第一行像元从左到有一一对应，亮度矩阵中的第二行值，与图象中的第二行像元从左到有一一对应，依次类推，直到最后一行。
在IDL中打开band1的亮度矩阵并且显示在IDL绘图环境中的时候，要正确地显示图象（见图2），必须要确保亮度矩阵与图象坐标系统的正确关系，如亮度矩阵中的左下三角（或右上三角）一定要严格对应图象像元矩阵的左下三角（或右上三角），否则就会发生图象翻转、倒置（见图3），其实质就是将波段亮度矩阵进行了转置或行列倒置。
练习：
程序调用。调用band4的ASCII文件的代码如下：
file=filepath('band1_matrix.txt',subdir='examples\data')
openr,lun,file,/get_lun
data=lonarr(300,200) ;事先知道图象大小
readf,lun,data
print,size(data,/n_elements) ；data就是band4图象的亮度矩阵
free_lun,lun

3.波段运算（Band Math）
Band MathTM 功能允许你处理导致单个波段输出的复杂表达式。这些数学表达式也可以应用于一个多波段文件中的所有波段，providing “File Math”。
关于使用波段运算的更多信息，请参阅 ENVI Programmer’s Guide 第 29 页的 “Band Math Basics”。
1．可利用的波段运算功能（Available Band Math Functions）
Band Math 功能为用户提供一个灵活的图像处理工具，其中许多功能是无法在任何其它的图像处理系统中获得的。该功能的能力与 IDL 语言的能力直接相关。可用的函数包括但不仅限于 表 4-2 中列出的数学表达式。
表 4-2: 一些可用的波段运算函数。
Series and Scalar 数学 三角函数 其它波段运算选项
加 （+） 正弦 （sin（x）） 关系运算符 （EQ、NE、LE、LT、GE、GT）
减 （-） 余弦 （cos（x）） 逻辑运算符 （AND、OR、XOR、NOT）
乘 （*） 正切 （tan（x）） 类型转换函数（byte， fix， long， float， double， complex）
除 （/） 反正弦 （asin（x）） IDL 返回数组结果的函数
最小运算符 （<） 反余弦 （acos（x）） IDL 返回数组结果的程序
最大运算符 （>） 反正切 （atan（x）） User IDL 函数和程序
绝对值 （abs（x）） 双曲正弦 （sinh（x））
平方根 （sqrt（x）） 双曲余弦 （cosh（x））
指数 （^） 双曲正切 （tanh（x））
自然指数 （exp（x））
自然对数 （alog（x））
以10为底的对数 （alog10（x））
注意
一些有效的 IDL 表达式要求整个输入数组存在于内存中，它可以不必与 ENVI tiling 操作相兼容。
2．Band Math 对话框
（1）. 选择 Basic Tools > Band Math.
将出现 Band Math 对话框。假如运算结果是一个二维数组，它将接受任何有效的 IDL 数学表达式、函数或程序。
（2）. 在标签为 “Enter an expression:” 的文本框内，输入变量名（将被赋值到整个图像波段或可能应用到一个多波段文件中的每个波段） 和所需要的数学运算符。
变量名必须以字符 “b” 或 “B” 开头，后面跟着 5 个以内的数字字符。
实例：
若你想计算三个波段的平均值，则在文本框“Enter an expression:”内输入数学方程式：
（（float（b4）-float（b3））/（（float（b4）+ float（b3））
这时，变量b4、b3自动跳入”Previous band math expression”对话框中，可以输入到文本框中。该表达式中使用的三个变量，“b3” 是第一个变量，“b4” 是第二个变量，注意，在本例中，IDL 的浮点型函数用来防止计算时出现字节溢出错误。
（3）. 输入一个有效的表达式，点击 “OK”处理。
将出现 Variable/Band Name Pairings 对话框。请参见以下部分。
4．要重新使用、保存或取消任何以前应用的数学表达式：
（1）. 点击显示在 “Previous Expression:” 列表中的任何表达式，把它导入到 “Enter an expression:” 文本区中。
（2）. 一旦被导入，点击 “OK”，把该表达式应用到一组新的波段。
将出现 Variable/Band Name Pairings 对话框。请见下列的详细向导。
5． 要把表达式保存到一个输出文件，点击 “Save”，然后当出现 Enter Output Filename 对话框时，键入输出文件名。
为了保持一致，输出文件名应该指定扩展名为 .exp 。
（Tips:可以在记事本中先输入复杂的公式，再导入公式文件）。
6．要恢复原先保存的表达式，点击 “Restore”，然后选择适当的文件名。
该表达式将显示在 “Previous Expression:” 列表中。
要清除所有原先的表达式，点击 “Clear”。

图 4-13: Band Math 对话框。
Variable/Band Name Pairings 对话框
Variable/Band Name Pairings 对话框允许你从一个输入波段列表中，把波段赋值给输入在 “Enter an expression:” 文本框中的变量。
要把一个值赋给原先实例中的变量 “b1”：
1. 在标签为 “Variables used in expression:” 的文本框内，点击表达式 “B1”。
2. 在标签为 “可利用波段列表:” 的列表中，点击所需要的波段。
注意，一旦第一个波段被选择，只有那些相同空间大小的波段被显示在波段列表中。
3. 按照同种方法，为 “B2”、“B3” 等赋予一个值。
要把一个多波段图像赋值给一个或所有变量：
1. 点击 “Map variable to Input file”。
2. 使用标准的 ENVI 文件选择步骤，选择一个文件（这可视为 “File Math”）。
所选择的文件可以是波谱子集，但是若一个以上的文件被使用，它们必须有相同的波段数。
通过数学表达式修改的文件数学（file math），一个多波段输出图像产生。
3. 一旦所有变量被定义，标准的 ENVI 输出对话框显示在 Variable/Bands Pairings 对话框的底部。可以把一个波段运算结果作为新的变量输入下一个运算过程。
要选择一个空间子集：
1. 点击 “Spatial Subset”。
2. 将出现标准的 File Spatial Subset 对话框 （系统默认值被设置为处理整个空间场景。
要把结果输出到一个文件或内存，选择 “File” 或 “Memory” 切换按钮。
若选择输出到一个文件，键入一个输出文件名，或使用 “Choose” 按钮选择一个文件名，然后点击 “OK”。结果图像被显示在可利用波段列表中。
将用公式计算的NDVI和用菜单命令直接计算的NDVI进行比较。
4.运用 IDL 程序进行波段运算
由于 ENVI 为你提供对 IDL 性能的访问，你可以使用内置的 IDL 功能部件的能力、IDL 用户函数，或书写你自己的程序执行自定义的操作。这些函数的唯一要求是它们接受一个或多个图像阵列作为输入，并且输出一个单波段二维数组的计算结果。这些函数必须保存在 IDL 路径列表内的一个目录下，以便它们将自动编译。通过使用 ENVI 主菜单 System 下拉菜单下的 Compile Module 选项（见 ENVI Programmer’s Guide 第 23 页的 “Incorporating New Routines”），也可以对它们进行编译。以下是用户波段运算功能的一些简单的实例。
波段运算函数 1
下面的实例是一个非常简单的自定义波段运算函数，它把两个波段相加。下面的程序文本可以在一个文本编辑器中输入，并用文件名 user_bm1.pro 来保存：
实例：
FUNCTION user_bm1， b1， b2
RETURN， b1+b2
END
要从 Band Math “Enter an expression:” 文本框中调用该函数，使用语法：
user_bm1（b1，b2）
方法：保存以上代码——编译——在Band Math “Enter an expression:”中输入user_bm1（b1，b2）——指明b1,b2对应的波段，OK
波段运算函数 2
下面的实例是一个自定义的波段运算函数，它把一个变量的数据类型转换为字节型，并将数值倒置（inverts the values）。下面的程序文本可以在一个文本编辑器中输入，并用文件名 user_bm2.pro保存：
实例：
FUNCTION user_bm2， b1
lut = 255 - BINDGEN（256）
b1 = BYTSCL（b1）
b1 = lut（b1）
RETURN， b1
END
要从 Band Math “Enter an expression:” 文本框中调用该函数，使用语法：
user_bm2（b1）
波段运算函数 3
下面的实例是一个自定义的波段运算函数，当 b1非零时，它用变量 b2 的值代替变量 b1。这一函数对分类图象非常有用，它用于将另一幅图象的像元代替分类的像元。下面的程序文本可以在一个文本编辑器中输入，并用文件名 user_bm3.pro 来保存：
实例：
FUNCTION user_bm3， b1， b2
b1 = （b1 NE 0）*b2
RETURN， b1
END
要从 Band Math “Enter an expression:” 文本框中调用该函数，使用语法：
user_bm3（b1，b2）
通过。
波段运算函数 4
下面的实例是一个自定义的波段运算函数，它计算归一化差值植被指数（Normalized Difference Vegetation Index，NDVI），并把它缩放到字节数据范围。注意，“min” 和 “max” 关键字在函数中是必需的，以确保同样的最小和最大值被用于缩放一个tiled图像中所有的tiles。
对于变量 b1，应该使用一个 0.8 ?m 附近的红外图像波段，而对于变量 b2，应该使用一个 0.6 ?m 附近的“红”波段。下面的程序文本可以在一个文本编辑器中输入，并保存为 user_bm4.pro：
实例：
FUNCTION user_bm4， b1，b2
b1=bytscl （（float（b1） - b2） / （float（b1）+b2）， min=-1.0，max=1.0）
RETURN， b1
END
要从 Band Math “Enter an expression:” 文本框中调用该函数，使用语法：
user_bm4（b1，b2）
注意显示技巧：显示为RGB时，将r、g、b分别赋予band3、ndvi、band3，效果明显。

实例5：
作用：The following example creates a custom Band Math function to perform the following ratio and optionally check for divide-by-zero.
步骤：在记事本或IDL中编写如下代码——文件名为”mfband.pro”——保存在ENVI中的“SAVE_ADD”目录下.——重启ENVI——调入含至少两个波段的图象——输入变量b1,b2——在"Enter Expression" text box中输入表达式bm_ratio（b1, b2）
function bm_ratio, b1, b2, check=check
den = float(b1) - b2
if (keyword_set(check)) then ptr = where(den eq 0., count) $
else count = 0
if (count gt 0) then den[ptr] = 1.0
result = (float(b1) + b2) / den
if (count gt 0) then result[ptr] = 0.0
return, result
end