操作

基于RSGIS监测洪灾变化上机操作实例

基本原理：

① 大气校正

遥感图像在获取过程中，受到大气吸收与散射、传感器定标、地形等因素的影

响，且会随时间的不同而有所差异。利用多时相遥感图像的光谱信息检测地物变化

的重要前提是要消除不变地物的辐射值差异。

大气校正的目的是消除大气和光照等因素对地物反射的影响，大多数情况下，

大气校正是反演地物真实反射率的过程。

目前可以进行大气校正的模块有很多种，如最早的 MODTRAN 4+，6S (Second

Simulation of the Satellite Signal in the Solar Spectrum)，ACORN，ATREM，在

ERDAS IMAGINE 8.7上的模块 ATCOR，以及 ENVI上的模块 FLAASH（基于

MODTRAN）。

FLAASH可对LANDSAT，SPOT，AVHRR，ASTER，MODIS，MERIS，AATSR，

IRS等多光谱、高光谱数据、航空影像及自定义格式的高光谱影像进行快速大气校

正分析。下面的大气纠正步骤，都是基于FLAASH进行的。

② 辐射定标

当我们拿到一幅原始影像，先要进行辐射定标，目的是把图像上的DN（Digital

Number）值转为辐亮度或者是反射率。辐射定标的结果可以是表观辐亮度（L），

也可以是表观反射率（ρ）。

计算表观辐亮度（L）的公式为：

Radiance=（（Lmax-Lmin）/（Qcalmax-Qcalmin）*（Qcal-Qcalmin）+Lmin ①

其中：Radiance 是表观辐亮度，注意单位是W/m2·sr·μm；Qcal为像元DN

值（也就是影像数据本身）； Qcalmax为传感器处最大辐亮度值所对应的DN值，

一般为255；Qcalmin 为传感器处最大辐亮度值所对应的DN值，一般为0；Lmax

和Lmin是从参数表中查询，Lmin为光谱辐亮度的最小值，单位同L；Lmax为光

谱辐亮度的最大值，单位同L。

计算表观反射率（ρ）的公式为：

ρ =π*L*d

2

/（ESUN*cos（θ）） ②

其中：ρ为表观反射率；L为①式中计算出来的表观辐亮度；d为日地距离；

ESUN为大气层外的太阳辐射，也可以说是传感器接收处的太阳辐射；θ为太阳天

顶角（这个可以通过影像的元数据获取）。以上参数可以查询下表获得。

对于TM影像来说，当LANDSAT传感器获取影像以后（Level 0），会将其转

化为 32 位浮点型的绝对辐亮度。之后进一步处理，将绝对辐亮度变为8位的DN

值，这也就是我们购买后拿到的数据（Level 1）。

如果要将 L1的DN值转化为传感器处的辐亮度值（at-nsor spectral radiance），

对于TM是8bit来说，Qcalmax取255，Qcalmin取0，可以将公式①简化为下面这

个公式：

上面的这个公式也可以改为：

其中，

Grecale即为gain，Brescale即为offt。各个波段的，以及

和见下表。

需要注意的是，上述参数在2003年5月5日前后是不一致的，所以在操作时，

一定要搞清楚影像获取的时间。ETM+影像与TM影像的Lmin、Lmax、Grescale、

Brescale、d及ESUN的参数不一样，可以查询下列表格：

ETM+ Spectral Radiance Range

watts/(meter squared * ster * µm)

Before July 1, 2000 After July 1, 2000

Band

Number

LMIN LMAX LMIN LMAX LMIN LMAX LMIN

Low Gain High Gain Low Gain High Gain

LMA

X

191.6 1 -6.2 297.5 -6.2 194.3 -6.2 293.7 -6.2

196.5 2 -6.0 303.4 -6.0 202.4 -6.4 300.9 -6.4

152.9 3 -4.5 235.5 -4.5 158.6 -5.0 234.4 -5.0

157.4 4 -4.5 235.0 -4.5 157.5 -5.1 241.1 -5.1

31.06 5 -1.0 47.70 -1.0 31.76 -1.0 47.57 -1.0

12.65 6 0.0 17.04 3.2 12.65 0.0 17.04 3.2

10.80 7 -0.35 16.60 -0.35 10.932 -0.35 16.54 -0.35

158.3 8 -5.0 244.00 -5.0 158.40 -4.7 243.1 -4.7

ETM+ Solar Spectral Irradiances

Band watts/(meter squared * µm)

1 1969.000

2 1840.000

3 1551.000

4 1044.000

5 225.700

7 82.07

8 1368.000

Earth-Sun Distance in Astronomical Units

Julian Julian Julian Julian Julian

Day Day Day Day Day

1 .9832 74 .9945 152 1.0140 227 1.0128 305 .9925

Distance Distance Distance Distance Distance

15 .9836 91 .9993 166 1.0158 242 1.0092 319 .9892

32 .9853 106 1.0033 182 1.0167 258 1.0057 335 .9860

46 .9878 121 1.0076 196 1.0165 274 1.0011 349 .9843

60 .9909 135 1.0109 213 1.0149 288 .9972 365 .9833

③ 气溶胶反演

我们采用FLAASH模块中的暗目标法来反演气溶胶光学厚度。暗目标法是由

Kaufmaun提出，它是利用660nm和2100nm处的反射率估算气溶胶量，由于2100nm

波长比大部分气溶胶微粒的直径要大，故该波段受气溶胶影响可以忽略；在大量的

试验中，他发现2100nm的植被反射率与660nm植被反射率之间存在稳定的关系，

因此可以直接利用2100nm的植被反射率来获取660nm处的植被反射率。气溶胶的

影响会使得实际获得的植被反射率与理论反射率存在一定差异，FLAASH模块中

正是利用了这个差异来反演气溶胶的光学厚度值。

要在FLAASH模块中获取图像气溶胶含量，传感器必须具有660nm和 2100nm

附近的通道，这些通道主要是用于获取“黑暗像元”，条件为，

如果输入图像中还具有 800nm 和 420nm 附近的通道，可以用于消除阴影和水体，

条件为。

数据准备：

经过几何粗校正、没有经过大气校正的Landsat-TM5（30米分辨率）遥感影像

两景，分别位于2010-7-2与2010-9-4文件夹中。

其中有关影像的相关信息存放于后缀为_的文件中，TM影像1-5及7

波段的中心波长信息，存放于.txt文件中，用于大气校正。

中心波长

处理过程：

①大气校正：打开ENVI软件，点击主菜单栏的file下的open image file，打

开2010-9-4中的LT5152KHC00_B1—B5及B7影像，再选择主菜单栏

的Basic Tools---Preprocessing---Calibration Utilities---Landsat Calibration，对6个波

段进行辐射定标。

在Landsat Calibration Input File列表中，先选择LT5152KHC00_B1

第一波段影像进行处理，点击OK。

参考处理影像相关信息文件LT5152KHC00_修改ENVI

Landsat Calibration对话框参数，如下图所示，输出文件命名为band1_。按照

相同方法依次对其余5个波段进行辐射定标，处进行波段选择

时需要注意。

选择主菜单栏的File---Save File As---ENVI Standard，点击Import Files，导入

辐射定标后的band1_—band5_及band7_ 6个波段文件，再点击

Reorder Files对6个波段从1至7依次进行排序，最终保存为的20100904_ca文件

包含6个波段。

接下来进行FLAASH校正操作。选择主菜单栏的Spectral---FLAASH，出现

FLAASH Atmospheric Correction Model Input Parameters窗口，FLAASH 模块的操

作界面分为三块：最上部设定输入输出文件；中间设定传感器的参数；下部设定大

气参数。

首先设定输入输出文件。FLAASH 模块要求输入辐亮度图像，输出反射率图

像。之前我们进行了辐射定标，得到辐亮度图像，在这里要把 BSQ 格式的图像转

换为 BIL 或者BIP 格式的图像，选择主菜单栏的Basic Tools---Convert Data(BSQ，

BIL，BIP)。然后再 Input Radiance Image 中选择转换格式后的图像。

当输入图像后，需要导入已经准备好中心波长.txt文件，其中含有一列TM 每

个波段中心波长的信息。 然后需要选择Scale Factor，即原始辐亮度单位与ENVI 默

认辐亮度单位之间的比例。ENVI 默认的辐亮度单位是μW/cm·sr·nm，而之前

2

我们做辐射定标时单位是W/m·sr·μm，二者之间转换的比例是10，因此在下图

2

中选择Single scale factor，填写10.000。

在Output Reflectance File里面设定输出文件的文件名和位置。

设定传感器参数。首先是Scene Center Location，即遥感图像中心的坐标，以

及Flight Date、Flight Time GMT（the Greenwich Mean Time），后两者可以在TM 的

LT5152KHC00_文件中找到，填入即可。 遥感图像中心坐标可

以在TM的任一波段的头文件中获取到，在任一波段文件处右键点击Edit Header，

进入后点击Edit Attributes---Map Info，点击，可以选择DMS或DDEG中任一

种形式填写入FLAASH窗口。

在Sensor Type 菜单中选择Landsat TM5，此时Sensor altitude 自动填上为

705km。而Pixel Size 填为30m。

根据遥感影像研究区实际情况，填写Ground Elevation，因为巴基斯坦南部地

区平均海拔在20m左右，所以填写为0.02km。

最关键的为大气参数部分:

a) Atmospheric Model(大气模式): 共有Sub-Arctic Winter (SAW)，Mid-Latitude

Winter (MLW)，U.S. Standard (US)，Sub-Arctic Summer(SAS)，Mid-Latitude

Summer (MLS) 和Tropical (T)。根据经纬度和时间可以选定研究区的大气模式，

见下表，这里选Tropical (T)。（研究属于北纬28°）

b) Aerosol Model（气溶胶模式）：有Rural，Urban，Maritime和Tropospheric四种

选择。根据实际情况选择即可。关于此四种模式的解释见下图。

c) 这里选Maritime。

c) 当我们选择TM 时，可选的参数还有Aerosol Retrieval 和Initial Visibility。这两

个参数对最后的结果有相当重要的影响，因此最好能调查到当地的Initial Visibility。

此外，AERONET 在全世界各地有测定AOD（Atmospheric Optical Depth）的站点，

可以查询AOD 以后转换为消光系数，通过消光系数估算能见度，此步骤比较繁琐，

在此不予详述。如果采用Aerosol Retrieval 中的K-T算法计算Visibility，且能够计

算出结果的话，则采用K-T 算法的能见度。

d) 关于Aerosol Retrieval。如果选择了下拉菜单中的K-T method，那么需要在

Multispectral Settings 中设定参数，在Kaufman-Tanre Aerosol Retrieval/Assign

Default Values Bad on Retrieval Conditions 中选择Over-land Retrieval Standard

(660:2100nm)即可。根据不同的研究区可以设定不同的模式。

这些参数用于确定黑暗像元，用于气溶胶反演；

KT Upper Channel：建议选择 2100nm附近的通道；

KT Lower Channel：建议选择 660nm附近的通道；

Maximum Upper Channel Reflectance：建议设置为 0.1，即：。

Reflectance Ratio：为反射率比值，建议设置为0.45，即：

Cirrus Channel (optional)：确定云的通道，建议设置为1367—1383nm左右的通

道。

水汽、气溶胶和云反演通道设置波长范围如下图所示：

上面这些参数波长范围的选择可以参考中心波长.txt文件，设置完Apply即可。

FLAASH校正完后，生成校正后的影像及结果文档，可以看到可见度和平均

水量结果。

对比经过大气校正与校正前的影像，发现经过大气校正后的影像，植被、水体

与裸地波谱曲线更能真实地反映情况。

植

被

水体

裸地

对2010-9-4影像进行大气校正后，根据相同原理及步骤对2010-7-2影像进行

大气校正（叙述省略，可以参照截图）。

① 几何校正：因为两景影像都是从USGS网站下载的经过几何粗校正的，所

以考虑不用对其进行几何校正这一过程。

② 提取水体范围：运用下式提取巴基斯坦洪灾部分区域的水体面积。

MNDWI=（TM2- TM5）/（TM2+ TM5）

选择主菜单栏的Basic Tools---Band Math，点击Import Files，输入计算公式：

float((b1-b2))/float((b1+b2))

点击OK。其中b1选择TM2波段，b2选择TM5波段，命名保存至文件路径，

点击OK进行计算。

通过计算得到的影像中，水体为正值，在0-1之间，土壤、植被及建筑物等物

体为负值。

接下来继续使用Band Math功能对提取MNDWI后的影像进行二值化，提取出

水体范围。因为水体的值大于等于0，而其他地物的值小于0，所以可以使用下列

语句对水体赋值为1，而其他地物赋值为0：

(b1 ge 0)*1+(b1 lt 0)*0

其中ge 、lt分别表示“大于等于”和“小于”。

点击OK。其中b1选择计算好的MNDWI波段，命名保存至文件路径，点击

OK进行计算。

计算得到二值化后的影像。

按照相同原理和方法对2010-7-2影像进行MNDWI提取和二值化，步骤省略。

因为影像最左侧边缘有影像质量问题引起的干扰，所以需要对其进行裁剪，以

免造成误差。这里我们选用感兴趣区进行裁剪的方法。选择主菜单栏的Basic

Tools---Region Of Interest---ROI Tool。

我们可以较为粗略的画一个框覆盖住有误差的地区，画完后在框内点击右键结

束。

选择主菜单栏的Basic Tools---Subt Data via ROIs，用刚画好的方框对影像

进行裁剪。

③ 监测洪灾变化

在ENVI中，选择主菜单栏的File---Save File As---ERDAS IMAGINE，将二值

化后的两景影像另存为可以ARCGIS中打开的img格式。

然后在ARCGIS中打开，最初影像显现为灰色。右键单击影像名称，选择

properties，点击Symboloy便签栏，点击左侧Unique Values，使影像显示为二值

化分类方式，依照相同方法对另一景影像进行设置。

运用空间分析工具中Raster Calculator功能对两景影像进行减法运用，用

2010-9-4减去2010-7-2，得到巴基斯坦洪灾前后水体变化的范围及面积。

运算后将影像计算分为三类：-1代表水体减少的地区，0代表水体无变化的地

区，而1代表水体增多的地区。

右键点击Calculation，选择Open Attribute Table，打开其属性表，可以发现每

一类型的count数，因为影像的分辨率是30m，所以每一类型的面积area可以通过

count*栅格面积（30*30=900）计算得到。

所以巴基斯坦这一影像地区的水体变化面积情况为：

面积(km)

2

受灾地区（1） 无影响地区（0） 水体减少地区（-1）

8188.98 38868.65 467.59