ENVI下基於劈窗算法從MODIS數據中反演海表溫度


劈窗算法最初是為反演海面溫度開發的,具體地說是針對NOAA/AVHRR的4和5通道設計的,后來也被用來反演地表溫度,這種算法較成熟,精度也高。劈窗算法以地表熱輻射傳導方程為基礎,利用10~13μm 大氣窗口內,兩個相鄰熱紅外通道(一般為10.5~11.5μm、11.5~12.5μm)對大氣吸收作用的不同,通過兩個通道測量值的各種組合來剔除大氣的影響,進行大氣和地表比輻射率的修正。
表達式為:TS= T4+ A(T4- T5) + B
其中:TS為地表真實溫度,T4和 T5分別為AVHRR的4和5通道,A和B為常量。
AVHRR 的通道4(10.15~11.13μm) 和通道5 (11.15~12.15μm) 恰與MODIS 的第31 波段(10.178~ 11.128μm) 和32 波段 (11.177~ 12.127μm ) 的中心波長相對應, 可將MODIS 的31、32 波段資料, 用 於劈窗算法進行地表溫度計算。很多學者對這個算法進行了推演,得到很多新的算法,如覃志豪、毛克彪等人。本文就是使用其他學者推演的算法。
利用MODIS數據劈窗算法反演海表溫度技術流程如下圖:

  

              圖1 技術流程圖
注:(1)按照本流程反演出來的結果是SST。陸地上的值可以視為無效值,若要得到正確的陸表溫度,需要加入海陸分離的步驟,以及城鎮和自然表面的比輻射率計算。
(2)MODIS數據下載:nasa官網:http://modis.gsfc.nasa.gov/.
http://ladsweb.nascom.nasa.gov/data/search.html
下面詳細介紹處理流程。
第一步:數據打開
Open->Open as->EOS->MODIS,選擇MOD021KM.A2013185.0245.005.2013185094144.hdf文件,點擊OK打開。分為三個數據集:熱紅外數據(20-36波段),可見光到短波紅外的輻射率數據(1-19、26波段),可見光到短波紅外的反射率數據(1-19、26波段)。
第二步:輻射定標
由於ENVI默認的設置是對MODIS數據進行自動的輻射定標,所以打開的數據即是經過了輻射定標的數據。

圖2 打開modis02級數據
第三步:幾何校正
MODIS1b數據是hdf格式,自帶有經緯度坐標信息,可自動進行幾何校正,分別對熱紅外數據集和反射率數據集進行自動幾何校正。
反射率數據集合校正:
(1)打開工具/Geometric Correction/Georeference by Sensor/Georeference MODIS,選擇反射率數據集,點擊Spatial Subset,選擇大亞灣的大概位置,點擊Spectral Subet,選擇2和19波段。點擊OK。

圖3 選擇數據同時選擇空間子集和光譜子集
(2)在Georeference MODIS Parameter面板設置為UTM, WGS84,49帶。

圖4 Georeference MODIS Parameter面板
(3)點擊OK,在結果輸出面板,設置路徑和文件名輸出。
熱紅外數據集幾何校正:方法與反射率數據集校正方法一樣,選擇31和32兩個波段校正。
第四步:海表溫度反演
本文使用的是《高懋芳,覃志豪等,MODIS數據反演地表溫度的基本參數估計方法》中分裂窗算法模型進行海表溫度反演,旨在學習ENVI中的操作流程。
算法為:T s = A 0 + A 1*T 31 - A 2* T 32 ( 公式1)
其中:Ts 是 地 表 溫 度( K ) , T 31 和 T 32 分 別 是M OD I S 第 31 和 32 波段的亮度溫度; A 0, A 1 和 A 2 是分裂窗算法的參數,分別定義如下:
A 0 = [D 32( 1 - C31 - D 31) / ( D 32 C31 -D 31 C32) ] a31 - [ D 31( 1 - C32 - D 32) /( D 32 C31 - D 31 C32) ] a 32 (公式 2)
A 1 = 1 + D 31/ ( D32C31 - D 31 C32) + [ D 32( 1 -C31 - D 31) / ( D 32C31 - D31 C32) ] b31 ( 公式3)
A 2 = D 31/ ( D 32 C31 - D 31C32) + [ D 31( 1 - C32 - D 32) / ( D 32 C31 - D 31C32) ] b32 ( 公式4)
式中,a 31,b31,a 32 和 b 32 是常量, 根據 MODIS的波段特征確定, 在地表溫度 0 ~ 5 e 范圍內, 這些常量 分 別 可 取 a 31 = -64.60363 , b31 = 0.440817, a 32 = -68.72575, b32 = 0.473453。
上述公式的中間參數分別計算如下:
Ci = Ɛiτi (ɵ) ( 公式6)
Di = [ 1 -τi (ɵ)] [ 1 + ( 1 - Ɛ i) τi (ɵ)] ( 公式7)
式中: i 是指 MODIS的第31和32波段, 分別為 i= 31 或 32; τi (ɵ)是視角為 ɵ的大氣透過率;Ɛ i 是波段 i 的地表比輻射率。
由以上公式可以看出, 該算法要求衛星遙感器的31 和 32 波段數據來計算星上亮度溫度, 同時還要求已知大氣透過率和地表比輻射率, 才能進行地表溫度的反演。
下面是詳細操作步驟:
(1)大氣透過率計算
大氣透過率τi (ɵ)是計算地表溫度的基本參數, 通常是通過大氣水汽含量來估計。經過前人研究,可以用MODIS第 2 和 19 波段來反演大氣水分含量,然后再根據大氣水分含量與大氣透過率之間的關系來估計大氣透過率。對於MODIS圖像中的任何一個像元,其可能的大氣水分含量用下式估計

(公式8)
式中: ω是大氣水分含量( g*cm-2) , α和β常量,分別α= 0.02 和 β=0.6321;ρ19和ρ2分別是MODIS第19和2波段的地面反射率。
使用bandmath工具計算大氣水分含量:
表達式:((0.02-alog(b19/b2))/0.6321)^2
B19:第19波段反射率
B2:第2波段反射率
大氣透過率的計算中,水汽是最主要的考慮因素,毛克彪等將MODTRAN等大氣模型模擬出來的兩者的關系,應用到MODIS數據中,提高了地表溫度反演的精度和實時性,本文采用模擬效果較好的指數關系模擬方程,擬合度達到了0.99以上,公式為:
τ31= 2.89798-1.88366*exp{-[ω/(-21.22704)]} (公式9)
τ32= -3.59289+4.60414*exp{-[ω/(-32.70639)]} (公式10)
式中,ω是水汽含量。
使用bandmath工具計算大氣透過率:
31波段大氣透過率=2.89798-1.88366*exp(b1/21.22704)
32波段大氣透過率=-3.59289+4.60414*exp(b1/32.70639)
B1:大氣水分含量。
(2)地表比輻射率的估算
地表比輻射率主要取決於地表的物質結構,對MODIS來說,大致分為水面、城鎮和自然表面。對於反演來說,利用混合像元分解的方法,根據植被覆蓋率來計算自然表面和城鎮的比輻射率,水體的可以用常量:
Ɛ 31水體=0.996,Ɛ 32水體=0.992。
有了這些參數,我們就可以計算C31、C32、D31、D32中間參數,BandMath表示式分別為:
C31=0.996*b31
C32=0.992*b32
D31=(1-b31)*(1+(1-0.996)*b31)
D32=(1-b32)*(1+(1-0.992)*b32)
其中 B31:31波段大氣透過率
B32:32波段大氣透過率
(3)亮度溫度的計算
將圖像DN值定標維熱輻射強度之后, 可用Planck函數求解出星上亮度溫度, 計算公式如下:
T i = K i 2 / l n ( 1 + K i 1 /I i )
式中, K i 1和 K i 2 是常量,對於第 i = 31 波段, 分別為 K 31 , 1 = 729 .541636 W•m-2•sr-1•um-1 ,
K 31 , 2 = 1304.413871K ; 對於第 i = 32 波段, 為 K 32 , 1 = 474 . 6 84780 W•m-2•sr-1•um-1,, K 31 , 2 = 1196 . 978785 K。
使用bandmath工具計算31和32的亮溫。
31波段亮溫=1304.413871/alog(1+729.541636/b31)
B31: 31波段輻射亮度值
32波段亮溫=1196.978785 /alog(1+474.684780/b32)
B32:32波段輻射亮度值
(4)A 0, A 1 和 A 2參數計算
接下來我們計算A 0, A 1 和 A 2參數,Bandmath表達式分別為:
A0=b4*(1-b1-b3)/(b4*b1-b3*b2)*(-64.60363)-b3*(1-b2-b4)/(b4*b1-b3*b2)*(-68.72575)
A1=1+b3/( b4*b1-b3*b2)+b4*(1-b1-b3)/( b4*b1-b3*b2)* 0.44081
A2=b3/(b4*b1-b3*b2)+b3*(1-b2-b4)/(b4*b1-b3*b2)* 0.473453
其中,b1:C31
B2:C32
B3:D31
B4:D32
(5)溫度計算
把這些參數帶入公式1中計算溫度值,Bandmath表達式為:
T s=b0+b1*b31-b2*b32-273
其中:b0:A0參數
B1:A1參數
B2:A2參數
B31:B31亮溫值
B32:B32亮溫值
如下為得到的最終地表溫度反演結果,單位為攝氏度。通過/Statistics/Compute Statistics工具統計看到最初的結果會有一些數量很少的異常值,幾十個像素,這些異常值大多是處於影像的邊緣。如下圖為統計的結果,小於-40的像元有57個,大於32度的34個像元。可以將這些像元處理,如用以下bandmath處理:-40>b1<32。

圖5 初始結果統計

            圖6 最終反演結果


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM