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

1、單窗算法反演地表溫度教程1.1 算法原理1.1.1 單窗算法單窗算法（MW算法）是覃志豪于2001年提出的針對TM數(shù)據(jù)只有一個熱紅外波段的地面溫度反演算法。經(jīng)過眾多學(xué)者驗證，單窗算法具有很高的反演精度，且同樣適用于ETM+和landsat8數(shù)據(jù)。公式如下：Ts=（a6（1-C6-D6）也（146-口6）C6D6）Tsensor-D6Ta）/C6式中，LST為地表溫度（K）,Tsensor是傳感器上白亮度溫度（K）,Ta是大氣平均溫度（K）;a、b為參考系數(shù)，當(dāng)?shù)乇頊囟葹?-70C時，a=-67.355351b=0.458606C、D為中間變量，計算公式為：C=2J=（1-r）l+（l-f）r式

2、中，E為地表比輻射率，雷為地面到傳感器的大氣總透射率。因此單窗算法反演地表溫度的關(guān)鍵是計算得到亮度溫度Tsenso、地表比輻射率反大氣透射率下和大氣平均作用溫度Ta。1.1.2 參數(shù)計算 輻射亮溫計算利用Planck公式將圖像像元對應(yīng)傳感器輻射強(qiáng)度值轉(zhuǎn)換為對應(yīng)的亮度溫度值。公式如下式中，Tsenso為亮度溫度值；J影像預(yù)處理后得到的光譜輻射值,，2、單位為w/(m,sjNm),K1、K2為常量，可由數(shù)據(jù)頭文件獲取。計算圖像輻射亮溫之前，需采用輻射定標(biāo)參數(shù)將像元灰度值DN轉(zhuǎn)換為熱輻射強(qiáng)度值，公式如下：.AX-式中，Ml為增益參數(shù)，Al為偏移參數(shù)，該參數(shù)可直接在影像通文件數(shù)據(jù)中獲取，

3、且ENVI軟件中已經(jīng)集成，不需要自己在查找。 地表比輻射率計算根據(jù)覃志豪針對TM影像提出的混合像元分解法來確定區(qū)域地表福輻射率。對于城市區(qū)域，我們簡單的將其分為水體、自然表面和建筑表面三種，因此針對混合像元尺度上的地表比輻射率通過下式來估算：£=PR£+(1尸)R+/rV1'、1'7nimt式中，色為混合像元的地表比輻射率；Pv為植被覆蓋率；Rv為植被的溫度比率；Rm為建筑表面的溫度比率；£v表示植被法地表比輻射率，¥m表示建筑表面的地表比輻射率；*表示輻射校正項。根據(jù)覃志豪經(jīng)驗公式，v=0.986;%=0.972根據(jù)Ast

4、er提供的常用地物比輻射率光譜庫Landsat8數(shù)據(jù)特點(diǎn)，對于B10波段，fv=0.98672;后m=0.96767。對于B11波段，v=0.98990;%=0.977515。在地表相對較平整的區(qū)域，一般可取de=維地表起伏較大的區(qū)域,d£可根據(jù)植被的構(gòu)成比例估計。研究表明熱輻射相互作用在植被與裸土分別占一半時達(dá)到最大，所以可以根據(jù)如下經(jīng)驗公式來估計d£：當(dāng)Pvw0.5日Id&=0.0038Pv當(dāng)Pv>0.5時,d&=0.0038(Pv)當(dāng)Pv=0.5時,d£=0.0019RV和Rm的確定采用覃志豪等提出的經(jīng)驗公式計算：RV=0.9332+0

5、.0585*PVRm=0.9886+0.1287*PVPV為植被覆蓋率可用下式計算：PV=(NDVI-NDVIs)/(NDVIv-NDVIs)其中NDVI是歸一化植被指數(shù)，NDVIV、NDVIS分別是植被和裸土的NDVI值,由于沒有詳細(xì)的區(qū)域植被和土壤光譜資料，采用覃志豪經(jīng)驗值以NDVIV=0.70和NDVIS=0.05來進(jìn)行植被覆蓋度的近似估計。 大氣平均作用溫度的計算對于大氣平均溫度，覃志豪等也提供了一個經(jīng)驗公式，該公式是根據(jù)Modtran所提供的標(biāo)準(zhǔn)大氣推導(dǎo)出來的，如下表所示，其中T0為近地表溫度（K）。袤3T大氣平均作用溫度與地面附近氣溫的關(guān)系叫大氣的剖間類型兀與”的關(guān)系

6、美國1網(wǎng)6年平均大氣兀二二5.為96+0W的45To拉帶平均大氣丁17976£+0號1715T口中緯度夏季平均人氣Ta=160110+0.92621To中緯度冬季平均大氣T,-1$.2704+0J1118Z；般研究區(qū)在重慶的，選擇中緯度夏季平均大氣 大氣透射率的計算由于單窗算法對輸入的參數(shù)較敏感，因此需要較精確的大氣透射率參數(shù)來反演地表溫度。影響大氣透射率的因素較多，如氣壓、氣溫、氣溶膠含量、大氣水分含量、03、CO2、CO、NH4等對熱輻射傳導(dǎo)均有不同程度的作用，從而使地表的熱輻射在大氣中的傳導(dǎo)產(chǎn)生衰減。因此，準(zhǔn)確的大氣透射率求算比較復(fù)雜，需要較詳細(xì)的大氣剖面數(shù)據(jù)，但

7、在實際應(yīng)用中很難獲取詳細(xì)的大氣剖面數(shù)據(jù)。研究表明大氣水分含量對大氣透射率的變化影響較大，而其它因素由于很難在短時間內(nèi)發(fā)生較大變化而對大氣透射率的影響較小。因此，水分含量就成為大氣透射率估計的主要考慮因素。覃志豪等根據(jù)大氣中的水汽含量動態(tài)變化對大氣透射率的影響程度進(jìn)行了研究，并對大氣透射率進(jìn)行了估計，如表3.1。其中:w是大氣水汽含量。以重慶為研究區(qū)的研究一般選用是氣溫較高時的估算公式大氣透射率估計方程大氣剖面水分含里w/（gcm-2）高氣溫低氣溫0.4-1.61.6-3.00.4-1.61.6-3.0相關(guān)系大氣透射率估計方程數(shù)平方(R2)p6=0.97429O.08007w0.99611r6=

8、1.03141-0.11536w0.99827r6=0.9820O-0.09611w0.99463p6=1.0537914142w0.99899標(biāo)準(zhǔn)誤差(SEE)0.0023680.0025390.0033400.002375接下來只需要計算出大氣水分含量便可計算得到大氣透射率。學(xué)者楊景梅根據(jù)地面、高空氣象要輸資料，提出利用地面濕度計算整層大氣可降水量，由此計算得出整層大氣水汽含量的方法（楊景梅邱金恒，用地面濕度參量計算我國整層大氣可降水量及有效水汽含量的方法的研究，大氣科學(xué)，2002.01）公式如下：=q+q獷式中，為大氣水汽含量（g/cm2）,W為整層大氣可降水量（cm）；coC1為經(jīng)驗系

9、數(shù)，以重慶為研究對象，一般可選取成都的經(jīng)驗值來計算。C0=0.1274;C1=0.6878計算大氣可降水量公式如下：式中，W為整層大氣可降水量（cm）,e為地面水汽壓（hPa）,a0al為經(jīng)驗系數(shù)，可從以下公式計算得出（9之33。）（中3號）b03exp(-L?沂+2.74耳40.15)0.04exp(0.6H)-&+也0.17+&(<P>Ha)40.204(p<33.)0.051 (-25畫+0.25,0">20")丸-i-1-090020。)0.066""l(p-33)"+4.41式中，。為地理緯度（

10、°）,H為海拔高度。重慶市主城區(qū)緯度取29.5,平均高程取0.35km,則有。=0.19604e+0.04691其中：e是絕對水汽壓（hPa）。采用WMO在1996年公布的氣象儀器和觀測方法指南中公布的飽和水汽壓E公式。對于純水面用公式計算，適用于溫度-45-60CE=6.112exp17.6270/(243.12+T0)RH=e/E*100推算出:e=6.112exp17.62.1大氣麗而畚敝/*I同褊度F黑體福肺臣度計靠I刖表溫度“.1.2.1影像預(yù)處理(1)輻射定標(biāo)選擇RadiometricCorrection/RadiometricCalibration。在FileSelec

11、tion對話框中，選擇數(shù)據(jù)LC81230322013132LGN02_MTL_Thermal,單擊SpectralSubset選擇ThermalInfrared1(10.9),打開RadiometricCalibration面板。(2)大氣校正本教程選擇Flaash校正法。FLAASHAtmosphericCorrection,雙擊此工具，打開輻射定標(biāo)的數(shù)據(jù)，進(jìn)行相關(guān)的參數(shù)設(shè)置進(jìn)行大氣校正。注意：如果在多光譜數(shù)據(jù)輻射定標(biāo)時Scalefactor值忘記設(shè)置，可在本步驟中打開輻射定標(biāo)數(shù)時設(shè)置singlescalefaceor值為0.1,若已設(shè)置，則默認(rèn)值為1即T0/(243.12+T0)*RH式中

12、，RH為相對濕度，可以通過氣象觀測數(shù)據(jù)獲得，To為氣溫(C)。1.2地表溫度反演步驟軸射屋尿川-TQu大，校正-!NDviTrn/亮聯(lián)圖窿/*地表比附射串計端一XgtputDatsTyp?ScaleFcctcrL00一RedometricCi1i*，二krScalefactor不能改變，報錯。保持默認(rèn)1即可oJipnlvTIaASMSrilingsOutputFiLcmuc:II值DisplayresultOOKCancel1）InputRadianceImage打開輻射定標(biāo)結(jié)果數(shù)據(jù)；2）設(shè)置輸出反射率的路徑，由于定標(biāo)時候；3）設(shè)置輸出FLAASH校正文件的路徑，最優(yōu)狀態(tài)：路徑所在磁盤空間足夠

13、大；4）中心點(diǎn)經(jīng)緯度SceneCenterLocation自動獲取；5）選擇傳感器類型：Landsat-8OLI;其對應(yīng)的傳感器高度以及影像數(shù)據(jù)的分辨率自動讀??；6）設(shè)置研究區(qū)域的地面高程數(shù)據(jù)；7）影像生成時的飛行過境時間：在layermanager中的Lc8數(shù)據(jù)圖層右鍵選擇ViewMetadata,瀏覽time字段獲取成像時間；注：也可以從元文件“LC81230322013132LGN02_MTL.tXt中找到，具體名稱：DATE_ACQUIRED=2013-05-12;SCENE_CENTER_TIME=02:55:26.6336980Z;8）大氣模型參數(shù)選擇：Sub-ArcticSumm

14、er（根據(jù)成像時間和緯度信息選擇）；9）氣溶膠模型AerosolModel:Urban,氣溶膠反演方法AerosolRetrieval:2-band(K-T);10)其他參數(shù)按照默認(rèn)設(shè)置即可。11)多光譜參數(shù)設(shè)置中,K-T反演選擇默認(rèn)模式：Defaults->Over-LandRetrievalstandard(600:2100)波譜響應(yīng)函數(shù)：默認(rèn)指向ProgramFilesExelisENVI51classicfilt_funclandsat8_oli.sli把它重新指向：.ProgramFilesExelisENVI51resourcefilterfuncslandsat8_oli.

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

16、)植被覆蓋度計算第一步，利用多光譜遙感影像，計算研究區(qū)域歸一化植被指數(shù)NDVI,ENVI軟件中有工具集成，使用即可。第二步，利用波段運(yùn)算工具，根據(jù)覃志豪經(jīng)驗公式，計算研究區(qū)植被覆蓋度。(b1gt0.7)*1+(b1lt0.05)*0+(b1ge0.05andb1le0.7)*(b1-0.05)/(0.7-0.05)b1:NDVI(2)地表比輻射率計算第一步，計算de(b1le0.5)*0.0038*b1+(b1eq0.5)*0.0038*(1-b1)+(b1gt0.5)*0.0019b1:植被覆蓋度Pv第二步，計算RV與Rm,在波段運(yùn)算器中分別輸入下列式中：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:指的是de1.2.2 輻射亮溫計算K1K2需要根據(jù)不同影像查表得到：(132