基于單窗算法反演地表溫的ENVI操作教程

1、單窗算法反演地表溫度教程1.1 算法原理1.1.1 單窗算法單窗算法（MW算法）是覃志豪于2001年提出的針對TM數(shù)據(jù)只有一個熱紅外波段的地面溫度反演算法。經(jīng)過眾多學者驗證，單窗算法具有很高的反演精度，且同樣適用于ETM+和landsat 8數(shù)據(jù)。公式如下：式中，LST為地表溫度（K），Tsensor是傳感器上的亮度溫度（K），Ta是大氣平均溫度（K）；a、b為參考系數(shù)，當?shù)乇頊囟葹?-70時，a = -67.355351，b = 0.458606；C 、D為中間變量，計算公式為：式中，為地表比輻射率，為地面到傳感器的大氣總透射率。因此單窗算法反演地表溫度的關鍵是計算得到亮度溫度Tsenso、

2、地表比輻射率、大氣透射率和大氣平均作用溫度Ta。1.1.2 參數(shù)計算 輻射亮溫計算 利用Planck公式將圖像像元對應傳感器輻射強度值轉(zhuǎn)換為對應的亮度溫度值。公式如下式中，Tsenso為亮度溫度值；影像預處理后得到的光譜輻射值，單位為，K1 、K2為常量，可由數(shù)據(jù)頭文件獲取。計算圖像輻射亮溫之前，需采用輻射定標參數(shù)將像元灰度值DN轉(zhuǎn)換為熱輻射強度值，公式如下：式中，ML為增益參數(shù)，AL為偏移參數(shù)，該參數(shù)可直接在影像通文件數(shù)據(jù)中獲取，且ENVI軟件中已經(jīng)集成，不需要自己在查找。 地表比輻射率計算 根據(jù)覃志豪針對TM影像提出的混合像元分解法來確定區(qū)域地表福輻射率。對于城

3、市區(qū)域，我們簡單的將其分為水體、自然表面和建筑表面三種，因此針對混合像元尺度上的地表比輻射率通過下式來估算： 式中，為混合像元的地表比輻射率；PV為植被覆蓋率；RV為植被的溫度比率；RM為建筑表面的溫度比率；V表示植被法地表比輻射率，m表示建筑表面的地表比輻射率；d表示輻射校正項。 根據(jù)覃志豪經(jīng)驗公式，V=0.986；m=0.972。 根據(jù)Aster提供的常用地物比輻射率光譜庫Landsat 8數(shù)據(jù)特點，對于B10波段，V=0.98672；m=0.96767。對于B11波段，V=0.98990；m=0.977515。 在地表相對較平整的區(qū)域,一般可取d=0;在地表起伏較大的區(qū)域,d可根據(jù)植被的

4、構(gòu)成比例估計。研究表明熱輻射相互作用在植被與裸土分別占一半時達到最大,所以可以根據(jù)如下經(jīng)驗公式來估計d：當Pv0.5時,d=0.003 8Pv當Pv>0.5時,d=0.003 8 (1-Pv)當Pv=0.5時, d=0.0019RV和Rm的確定采用覃志豪等提出的經(jīng)驗公式計算:RV=0.9332+0.0585*PVRm=0.9886+0.1287*PVPV為植被覆蓋率可用下式計算:PV=(NDVI-NDVIS)/(NDVIV-NDVIS) 其中NDVI是歸一化植被指數(shù),NDVIV、NDVIS分別是植被和裸土的NDVI值,由于沒有詳細的區(qū)域植被和土壤光譜資料,采用覃志豪經(jīng)驗值以NDVIV=0

5、.70和NDVIS=0.05來進行植被覆蓋度的近似估計。 大氣平均作用溫度的計算對于大氣平均溫度，覃志豪等也提供了一個經(jīng)驗公式，該公式是根據(jù)Modtran所提供的標準大氣推導出來的，如下表所示，其中T0為近地表溫度(K)。一般研究區(qū)在重慶的，選擇中緯度夏季平均大氣。 大氣透射率的計算由于單窗算法對輸入的參數(shù)較敏感，因此需要較精確的大氣透射率參數(shù)來反演地表溫度。影響大氣透射率的因素較多,如氣壓、氣溫、氣溶膠含量、大氣水分含量、O3、CO2、CO、NH4等對熱輻射傳導均有不同程度的作用,從而使地表的熱輻射在大氣中的傳導產(chǎn)生衰減。因此,準確的大氣透射率求算比較復雜,需要較

6、詳細的大氣剖面數(shù)據(jù)，但在實際應用中很難獲取詳細的大氣剖面數(shù)據(jù)。研究表明大氣水分含量對大氣透射率的變化影響較大,而其它因素由于很難在短時間內(nèi)發(fā)生較大變化而對大氣透射率的影響較小。因此,水分含量就成為大氣透射率估計的主要考慮因素。覃志豪等根據(jù)大氣中的水汽含量動態(tài)變化對大氣透射率的影響程度進行了研究，并對大氣透射率進行了估計，如表3.1。其中:w是大氣水汽含量。以重慶為研究區(qū)的研究一般選用是氣溫較高時的估算公式。大氣透射率估計方程大氣剖面水分含量w/(g·cm-2)大氣透射率估計方程相關系數(shù)平方(R2)標準誤差(SEE)高氣溫0.996110.0023680.998270.002539低氣

7、溫0.994630.0033400.998990.002375接下來只需要計算出大氣水分含量便可計算得到大氣透射率。學者楊景梅根據(jù)地面、高空氣象要輸資料，提出利用地面濕度計算整層大氣可降水量，由此計算得出整層大氣水汽含量的方法（楊景梅 邱金恒，用地面濕度參量計算我國整層大氣可降水量及有效水汽含量的方法的研究，大氣科學，2002.01）公式如下：式中，為大氣水汽含量（g/cm2），W為整層大氣可降水量（cm）；C0 C1為經(jīng)驗系數(shù)，以重慶為研究對象，一般可選取成都的經(jīng)驗值來計算。C0= 0.1274；C1=0.6878.計算大氣可降水量公式如下：式中，W為整層大氣可降水量（cm），e為地面水汽壓

8、（hPa），a0 a1為經(jīng)驗系數(shù)，可從以下公式計算得出。式中，為地理緯度（°），H為海拔高度。重慶市主城區(qū)緯度取29.5，平均高程取0.35 km，則有。=0.19604e+0.04691其中：e是絕對水汽壓（hPa）。采用WMO在1996年公布的氣象儀器和觀測方法指南中公布的飽和水汽壓E公式。對于純水面用公式計算，適用于溫度-4560。 E=6.112exp17.62*T0/(243.12+T0) RH=e/E*100推算出： e=6.112exp17.62*T0/(243.12+T0)*RH式中，RH為相對濕度，可以通過氣象觀測數(shù)據(jù)獲得，T0為氣溫()。1.2 地表溫度反演步驟1

9、.2.1 影像預處理（1）輻射定標選擇Radiometric Correction/Radiometric Calibration。在File Selection對話框中，選擇數(shù)據(jù)，單擊Spectral Subset選擇Thermal Infrared1（10.9），打開Radiometric Calibration面板。Scale factor 不能改變，否則后續(xù)計算會報錯。保持默認1即可。（2）大氣校正本教程選擇Flaash 校正法。FLAASH Atmospheric Correction，雙擊此工具，打開輻射定標的數(shù)據(jù)，進行相關的參數(shù)設置進行大氣校正。注意：如果在多光譜數(shù)據(jù)輻射定標時S

10、cale factor值忘記設置，可在本步驟中打開輻射定標數(shù)時設置single scale faceor 值為0.1，若已設置，則默認值為1即可。1)   Input Radiance Image：打開輻射定標結(jié)果數(shù)據(jù)；2)  設置輸出反射率的路徑，由于定標時候；3)  設置輸出FLAASH校正文件的路徑，最優(yōu)狀態(tài)：路徑所在磁盤空間足夠大；4)  中心點經(jīng)緯度Scene Center Location：自動獲??；5)   選擇傳感器類型：Landsat-8 OLI；其對應的傳

11、感器高度以及影像數(shù)據(jù)的分辨率自動讀??；6)  設置研究區(qū)域的地面高程數(shù)據(jù)；7)   影像生成時的飛行過境時間：在layer manager中的Lc8數(shù)據(jù)圖層右鍵選擇View Metadata，瀏覽time字段獲取成像時間；注：也可以從元文件中找到，具體名稱：DATE_ACQUIRED = 2013-05-12；SCENE_CENTER_TIME = 02:55:26.6336980Z；8)  大氣模型參數(shù)選擇：Sub-Arctic Summer（根據(jù)成像時間和緯度信息選擇）；9)   氣溶膠模型

12、Aerosol Model：Urban，氣溶膠反演方法Aerosol Retrieval：2-band（K-T）；10)  其他參數(shù)按照默認設置即可。11) 多光譜參數(shù)設置中，K-T反演選擇默認模式：Defaults->Over-Land Retrieval standard（600:2100）波譜響應函數(shù)：默認指向.Program FilesExelisENVI51classicfilt_funclandsat8_oli.sli把它重新指向：.Program FilesExelisENVI51resourcefilterfuncslandsat8_oli.sli 注

13、：這是因為ENVI5.1版本的一個小bug，即Classic中的L8的波譜響應函數(shù)不正確，另外一個一勞永逸的方法是：將“ExelisENVI51resourcefilterfuncs”中的Landsat8_oli.sli 和 Landsat8_oli.hdr兩個文件拷貝覆蓋：“.ENVI51classicfilt_func” 中的兩個文件。否則SWIR1波段大氣校正后的結(jié)果全為0。見下圖12) 高級參數(shù)設置：根據(jù)內(nèi)存大小設置Tile Size（Mb）：100(8g物理內(nèi)存)，其他參數(shù)默認即可，詳細見下圖經(jīng)過FLAASH校正的影像基本去除了空氣中水汽顆粒等因子的影響，植被的波譜曲線趨于正常。1.

14、2.1.1 地表比輻射率計算（1）植被覆蓋度計算第一步，利用多光譜遙感影像，計算研究區(qū)域歸一化植被指數(shù)NDVI，ENVI軟件中有工具集成，使用即可。 第二步，利用波段運算工具，根據(jù)覃志豪經(jīng)驗公式，計算研究區(qū)植被覆蓋度。(b1 gt 0.7)*1+(b1 lt 0.05)*0+(b1 ge 0.05 and b1 le 0.7)*(b1-0.05)/(0.7-0.05)b1：NDVI（2）地表比輻射率計算 第一步，計算d(b1 le 0.5)*0.0038*b1+(b1 eq 0.5)*0.0038*(1-b1)+(b1 gt 0.5)*0.0019b1：植被覆蓋度Pv 第二步，計算RV與Rm，

15、在波段運算器中分別輸入下列式中：RV : 0.9332+0.0585*b1Rm : 0.9886+0.1287*b1b1: 植被覆蓋度Pv 第三步，計算地表比輻射率利用B10波段計算地表比輻射率: b1*b2*0.98672+(1-b1)*b3*0.96767+b4利用B11波段計算地表比輻射率: b1*b2*0.9899+(1-b1)*b3*0.97515b1: 植被覆蓋度Pv; b2:植被表面溫度比率RV; b3:建筑表面溫度比率Rm;b4：指的是d1.2.2 輻射亮溫計算K1 K2需要根據(jù)不同影像查表得到：（1321.0789）/alog（774.8853/b1 +1）b1：波段10 或