操作

基于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