本文介紹基於ENVI軟體，實現對Landsat 7遙感影像加以大氣校正方法的地表溫度反演操作。

1 影象前期處理與本文理論部分

更新：基於GEE的地表溫度Landsat反演可以看文章單窗演演算法的地表溫度反演：谷歌地球引擎GEE程式碼，自動批次操作，處理更快。

本文操作部分將直接由植被覆蓋度計算展開；而對於一個完整的地表溫度反演計算過程，在求解植被覆蓋度這一步驟之前仍有很多資料準備、預處理等工作。為了更好理解整個實驗過程，將我們未進行的步驟梳理如下。其中，具體的前期操作方法大家可以參考文章ENVI實現QUAC、簡化黑暗像元、FLAASH方法的遙感影像大氣校正。

首先獲取需要的遙感影象資料，並對其進行包括上述文章內容在內的預處理步驟——資料匯入、輻射定標、幾何校正、大氣校正、影象拼接與裁剪等。

其中，輻射定標需要分兩步驟，即對可見光波段資料（如上述文章中的1、2、3、4、5、7六個波段）與熱紅外波段資料（如上述文章中的6波段）分別進行輻射定標。

其次，需要計算NDVI（即Normalized Difference Vegetation Index，歸一化植被指數，而非植被覆蓋度）。NDVI是指一幅遙感影像中，近紅外波段的反射值與紅光波段的反射值之差比上這二者之和；其可以用來檢測植被生長狀態、植被覆蓋度，還可以消除部分輻射誤差等。NDVI的具體取值範圍限制在-1到1之間，其負值表示地面覆蓋為雲、水、雪等，即對可見光具有高反射；0值表示地面覆蓋有岩石或裸土等，從而使得NIR（Near Infrared，近紅外波段）和R（Red，紅光波段）近似相等；其正值則表示地面有植被覆蓋，且植被覆蓋度越高，其數值越高。目前，在一些網站（如NASA官方網站）具有NDVI成品資料，可供我們直接下載、利用；而通過初始遙感影像中的近紅外波段資料和紅光波段資料，我們可以直接利用前述定義公式，即

對其加以計算。計算NDVI時需要注意，所選用的遙感影像不能具有過多的雲干擾。

再次，需要視實際情況對計算得出的NDVI結果影象進行重取樣。這是由於，在本文中需要多次利用「Band Math」工具對影象資料進行計算，而這一工具要求輸入的資料在影象解析度（即像元大小）及行、列數等方面完全一致。同時，我們即將使用的熱紅外資料（即Landsat ETM+第六波段資料），其解析度為60米；而計算得到的NDVI資料影象解析度為30米。因此，我們需要對解析度精度更高的NDVI資料影象加以重取樣處理，使得二者解析度一致。

重取樣功能可以通過ENVI軟體中選擇「Basic Tools」→「Resize Data (Spatial/Spectral)」加以實現。

目前主要的地表溫度單波段反演演演算法包括大氣校正法（又名輻射傳輸方程法，Radiative Transfer Equation，RTE）、單通道演演算法和單窗演演算法。本文我們使用大氣校正方法。大氣校正法的基本原理為：估計得到大氣對地表熱輻射的影響，然後將這一部分大氣影響由衛星感測器所接收到的熱輻射總量中減去，得到地表熱輻射強度；最後將地表熱輻射強度轉化為對應的地表溫度即可。

衛星感測器接收到的熱紅外輻射亮度值Lλ由三部分組成，其分別為：大氣向上輻射亮度L↑、地表物體的真實輻射亮度在經過大氣層後到達衛星感測器的能量、大氣向下輻射亮度在經過地面反射後的能量。因此，結合上述這一理論過程，可以用輻射傳輸方程來表示衛星感測器接收到的熱紅外輻射亮度值：

其中，ε為地表比輻射率，T_s為地表真實溫度，B(T_s )為地表在T_s這一真實溫度下的黑體熱輻射亮度，τ為大氣的熱紅外波段透過率。地表比輻射率ε又稱為發射率，而根據基爾霍夫定律，發射率與吸收率相等，則(1-ε)可以表示反射率。因此，[εB(T_s )τ]即為地表物體的真實輻射亮度在經過大氣層後到達衛星感測器的能量，而[(1-ε)L↓τ]則表示大氣向下輻射亮度在經過地面反射後的能量。

物體的比輻射率是物體向外輻射電磁波的能力表徵，是指在同一溫度下地表發射的輻射量與一黑體發射的輻射量的比值。其不僅依賴於地表物體的組成，而且與物體的表面狀態（如表面粗糙度）及物理性質（如介電常數、含水量等）有關，並隨著所測定的波長和觀測角度等因素變化。對地表比輻射率的精確定量測量難度較大，因此本文依據經驗法對地表比輻射率加以估計：

其中，F_C為植被覆蓋度。由此可知，我們需要同時計算出地表植被覆蓋度，用以確定地表比輻射率計算公式並參與計算。本文我們採取混合像元分解法求解植被覆蓋度。同上述地表比輻射率計算公式較為一致，我們依然將地表分為水體、植被與建築三個部分；其中，依據NDVI數值對這一類別加以具體區分。NDVI小於0時，認為地物為水體，植被覆蓋度為0；NDVI大於0.7時，認為地物為植被，植被覆蓋度為1；NDVI取值在[0,0.7]時，認為地物處於水體與植被之間，植被覆蓋度依據公式計算。針對不同的地物，計算植被覆蓋度：

最後，完成上述全部計算並依據輻射傳輸方程求得B(T_s )後，可以依據普朗克公式反函數求出地表真實溫度。公式為：

所得地表真實溫度單位為開爾文（K），我們需要將其轉換為常見的攝氏度。溫度單位的轉換可以在得到地表真實溫度影象後單獨計算，亦可以直接在上式中計算。

2 實際操作

2.1 植被覆蓋度計算

依據上述分析，首先我們需要藉助60米解析度的研究區NDVI資料影象，計算研究區的植被覆蓋度。

（1） 開啟ENVI Classic 5.3（64-bit）軟體，選擇「File」→「Open Image File」，在彈出的檔案選擇視窗中選擇「TM-NDVI-60m.img」檔案；點選「開啟」。

（2） 選擇「Basic Tools」→「Band Math」，在彈出的公式建立視窗中輸入本次實驗的第一個公式，即植被覆蓋度公式。輸入公式完成後，點選下方「Add to List」按鈕，即可將公式存入待選擇區內。為了減少後期不必要的工作量，可以每次編輯完成一個公式後點選「Save」按鈕，以將待選擇區內的公式儲存。

（3） 儲存公式完成後，點選左下角「OK」按鈕，即可開始公式的計算。在彈出的公式變數檔案選擇視窗中，將這一公式的變數「B1」選擇為我們剛剛新增的檔案「TM-NDVI-60m.img」。隨後，設定輸出檔案地址等資訊。

（4） 設定完成後，點選左下角「OK」按鈕，即可開始公式的執行。執行結束，將所得到的研究區植被覆蓋度結果影象匯入ENVI軟體中，顯示如下。

2.2 地表比輻射率計算

如前所述，我們完成植被覆蓋度的計算其實是為地表比輻射率計算做準備。同樣是將地物分為水體、城鎮與自然表面三個類別，依據經驗法，運用與植被覆蓋度形式類似的分段計算公式，為不同類別地物賦予不同的地表比輻射率計算公式。

考慮到在計算植被覆蓋度時使用的公式中已包含「B1」這個變數波段名稱，為了避免不同公式之間的變數相互混淆，因此在一開始我在每一個公式中都使用了不同的變數名稱——如在此處地表比輻射率計算公式中使用「B2」「B3」等。但通過後期的實驗過程發現，其實每一次公式的變數名稱即使一致也不會對實驗造成明顯的影響。

（1） 選擇「Basic Tools」→「Band Math」，在彈出的公式建立視窗中輸入本次實驗的第二個公式，即地表比輻射率公式。輸入公式完成後，點選下方「Add to List」按鈕，即可將公式存入待選擇區內。為了減少後期不必要的工作量，可以在編輯完成這一公式後點選「Save」按鈕，並將最近一次儲存的公式檔案覆蓋，以將待選擇區內的兩條公式統一儲存。

（2） 儲存公式完成後，點選左下角「OK」按鈕，即可開始公式的計算。在彈出的公式變數檔案選擇視窗中，將這一公式的變數「B2」依然選擇為我們最開始新增的影象檔案「TM-NDVI-60m.img」，並將這一公式的變數「B3」選擇為通過上述步驟獲得的植被覆蓋度結果影象檔案。隨後，設定輸出檔案地址等資訊。

（3） 設定完成後，點選左下角「OK」按鈕，即可開始公式的執行。執行結束，將所得到的研究區地表比輻射率結果影象匯入ENVI軟體中，顯示如下。

2.3 相同溫度下黑體輻射亮度值計算

正本文第一部分分析所得，已知研究區域地表比輻射率與熱紅外波段亮度，我們便可以計算相同溫度下黑體輻射亮度值。

（1） 在ENVI Classic 5.3（64-bit）軟體中選擇「File」→「Open Image File」，在彈出的檔案選擇視窗中選擇「TM6-rad-subset-jz-xiangfan.img」檔案；點選「開啟」。

（2） 選擇「Basic Tools」→「Band Math」，在彈出的公式建立視窗中輸入本次實驗的第三個公式，即黑體輻射亮度值公式。輸入公式完成後，點選下方「Add to List」按鈕，即可將公式存入待選擇區內。為了減少後期不必要的工作量，可以在編輯完成這一公式後點選「Save」按鈕，並將最近一次儲存的公式檔案覆蓋，以將待選擇區內的三條公式統一儲存。

（3） 儲存公式完成後，點選左下角「OK」按鈕，即可開始公式的計算。在彈出的公式變數檔案選擇視窗中，將這一公式的變數「B4」選擇為我們剛剛新增的影象檔案「TM6-rad-subset-jz-xiangfan.img」，並將這一公式的變數「B5」選擇為通過上述步驟獲得的地表比輻射率結果影象檔案。隨後，設定輸出檔案地址等資訊。

（4） 設定完成後，點選左下角「OK」按鈕，即可開始公式的執行。執行結束，將所得到的研究區相同溫度下黑體熱紅外波段輻射亮度值計算結果影象匯入ENVI軟體中，顯示如下。

2.4 地表真實溫度計算

由前述分析可知，通過上述步驟得到的黑體熱紅外波段輻射亮度並不是地表實際溫度，我們依然需要通過普朗克公式反函數實現二者之間的轉換關係。
與此同時，這一步驟得到的地表實際溫度結果單位為開爾文（K），並不是我們平日裡經常使用的攝氏度（℃）。因此，我們還需要實現溫度單位的轉換。

由開爾文溫度轉換為攝氏度只需要在原溫度基礎之上減去273.15即可，較為簡單，沒有必要單獨轉換。因此我選擇直接在這一步驟將溫度單位的轉換完成。

（1） 選擇「Basic Tools」→「Band Math」，在彈出的公式建立視窗中輸入本次實驗的第四個公式，即地表真實溫度公式。輸入公式完成後，點選下方「Add to List」按鈕，即可將公式存入待選擇區內。為了減少後期不必要的工作量，可以在編輯完成這一公式後點選「Save」按鈕，並將最近一次儲存的公式檔案覆蓋，以將待選擇區內的四條公式統一儲存。

（2） 儲存公式完成後，點選左下角「OK」按鈕，即可開始公式的計算。在彈出的公式變數檔案選擇視窗中，將這一公式的變數「B6」選擇為我們剛剛計算獲得的黑體熱紅外波段輻射亮度結果影象檔案，隨後，設定輸出檔案地址等資訊。

（3） 設定完成後，點選左下角「OK」按鈕，即可開始公式的執行。執行結束，將所得到的研究區地表真實溫度計算結果影象匯入ENVI軟體中，顯示如下。

（4） 在得到的結果影象任意位置處右鍵，選擇「Cursor Location/Value」或「Quick Stats」選項，均可以檢視影象各個像元的畫素資訊。在這裡，我使用「Cursor Location/Value」檢視影象像元資訊，得到如下所示的結果。

（5） 可以看到，其中像元畫素值（即「Data」值）出現了「302.902771」這一數值。而這一步驟所得到的結果為單位為攝氏度的地表實際溫度數值，不可能出現三百多的資料。因此，說明很可能在前面操作部分出現錯誤。

（6） 通過返回檢查發現，在其中計算地表比輻射率時，所輸入的公式出現了輸入錯誤。將這一錯誤糾正後，重新生成地表實際溫度影象，並對其進行統計檢查。

（7） 可以看到，經過修改後的地表真實溫度影象資料符合實際情況，可以認為錯誤已被排除。

2.5 影象匯出

由於在實驗後期需要製作專題地圖，以將溫度分為不同等級並比較不同地物的溫度特性，因此需要將ENVI中得到的地表實際溫度結果影象另存為「.tif」格式，從而方便利用ArcGIS系列軟體對其加以進一步分級、美化與出圖處理。

（1） 首先，我在ENVI影象視窗中進行影象的另存為操作。選擇「File」→「Save Image As」→「Image File」，在彈出的儲存設定視窗中選擇檔案格式為「TIFF/GeoTIFF」，並設定好結果影象檔案儲存路徑、儲存檔名等。

（2） 得到儲存的結果影象檔案後，將其新增進入AcrMap 10.2軟體中。此時發現，得到的圖層檔案在畫素上已被拉伸至0-255範圍。通過AcrMap的「識別」功能，發現像元畫素已全部被拉伸，即其原值均已成為上述0-255範圍的資料。因此，認為影象應該不會再通過相關操作恢復原有的溫度數值。隨後，嘗試用「設定柵格屬性」等工具加以調整，均以失敗告終。

（3） 通過查閱網路資源，看到有人指出這種情況是儲存了「影象」而不是「影像」；在這一提醒下，嘗試直接利用ENVI軟體左上角工具列中的「File」→「Save File As」→「TIFF/GeoTIFF」對影象加以儲存。在彈出的儲存設定視窗中選擇檔案，並設定好結果影象檔案儲存路徑、儲存檔名等。

（4） 將第二次儲存的「.tif」格式影象結果檔案匯入AcrMap 10.2軟體中，看到結果數值為正常狀態。

（5） 得到正常結果後，對兩種不同的儲存結果加以思考。本文所使用的軟體為ENVI Classic 5.3（64-bit）版本，而使用非Classic版本同樣發現這一問題。回顧在影象介面儲存圖片的過程，進一步發現，這樣的儲存應當是將影象格式直接設定成為了8bit整形儲存的格式，從而丟失其原有的畫素資料資訊與意義——即其只儲存了這幅影象的「外貌」，運用0-255之間並無實際意義的資料來表示影象的灰度；但其原有影象，即我們希望得到的影象，其資料、灰度等級應是溫度的表示。若將這張原有由不同溫度劃分出不同灰度等級的影象改變為由0-255之間數位劃分出不同灰度等級的影象，其自然失去了原有的溫度意義。

（6） 將正常影象匯入AcrMap 10.2軟體中後，同時發現其影像的黑邊無法去除——運用「識別」功能可以看到，黑邊部分並無資料，其原本就已均為「NoData」狀態。針對這一現象，可能是不同軟體在儲存、讀取資料檔案時出現的常見錯誤，無需在意，後期對溫度重新劃分等級後即可消失。

（7） 最終匯入影象及其資訊如下圖所示。其中，有兩個結果相近的圖層是因為在一開始地表實際溫度轉換過程中，我減去的數值為「273」而並非「273.15」。

3 專題地圖製作

匯入AcrMap 10.2軟體中的影象為像元畫素由22.3997至47.0569，若不對其重新劃分溫度區間處理，將會導致所成彩色影象各個顏色十分離散，不利於觀察、分析。

（1） 匯入後的影象自身不具有統計資訊，無法在「已分類」模組中對其各個數值劃分割區間。因此，需要首先使用工具對其資料資訊加以統計。點選「Data Management Tools.tbx」→「柵格」→「柵格屬性」→「計算統計資料」。

（2） 在彈出的設定視窗中選擇需要統計資料的圖層，其它專案不需要做調整。

（3） 點選「確定」，即可對圖層影象完成統計資料的計算。

（4） 在圖層名稱處右鍵，選擇「屬性」→「符號系統」→「已分類」，開啟「分類」按鈕，將分類方法選擇為「手動」，並依據要求設定四個「中斷值」——30，35，39和影象畫素的最大值47.0569。點選確定即可儲存間斷點設定。

（5） 回到圖層屬性視窗後，對四個範圍所表示的內容加以標註，以方便後期製圖。同時，針對溫度專題圖的特色，調節合適的配色方案。原本我準備使用紅色系作為配色，但發現出圖後整體區分效果並不是很好。因此選擇另一種對比相對較強的配色方案。

4 不同地物地表溫度對比分析

為了更加清晰地對比不同地物地表溫度的差異，藉助騰訊地圖中衛星地圖模組，以襄陽市漢江為對照，對比研究區域衛星地圖與專題地圖。得出結論如下：

（1） 水體溫度明顯低於周圍城市、鄉村等區域。這一特徵可謂溫度專題地圖中最為明顯的要素之一——由專題地圖可以清晰識別出一條低溫帶由襄陽市西部進入，由西北向東南方向延申；並在城市中部偏西南方向形成衝擊島，隨後繼續向南行進。此外，在衝擊島北部，由唐白河形成的低溫帶同樣較為明顯；其清晰可辨，甚至可以沿低溫帶追溯至該河流的發源處。城市東南處的秦咀水庫同樣如此。由水體的溫度色帶可以看出，在成像時刻，其溫度在30℃以下。由於水體的溫度和光照有密切關係，尤其是漢江、唐白河這一類地上河；而本次實驗所用影象的成像時間為上午10:30左右，當天日照時間還不長，因此水溫較低。

（2） 城鎮溫度明顯高於周圍區域。由溫度專題地圖可以看出，主要的高溫聚集區域均位於漢江、唐白河沿岸；結合衛星地圖可以看出，這些大多數均為襄陽市主要城市聚集中心，建築、人口密集，經濟發達。由城鎮溫度色帶可以看出，在成像時刻，其溫度在39℃以上。在另一方面，對於一些建築物分佈較為零散的郊區、農村等，溫度相對沒有城市中心那麼高；但其溫度普遍在35-39℃區間內，較之周圍區域溫度還是較高。由於城鎮人口多，汽車、生活燃燒等排放較多，加之工業生產頻繁，因此溫度較高。

（3） 農田、耕作用地等溫度較高。35℃以上區域主要分佈在漢江、唐白河沿岸的另一原因是這些耕作用地和城鎮一樣，多分佈於地勢平坦、土壤肥沃、灌溉方便的河流兩岸。由專題地圖可以看出，多數耕作用地的溫度在35-39℃這一區間；當然，同樣亦有不小面積的農田溫度在30-35℃區間。因此，可以認為耕作用地溫度處於35℃左右。其溫度較之城鎮低，是由於鄉村地廣人稀，人口密度顯著小於城市地區，因此產生的熱量較少；而其溫度較高於森林、山區則是因為地形平坦、人為耕作調節土溫等。

（4） 森林、山區等溫度明顯較低。在溫度專題地圖中，有三處主要的低溫聚集中心，分別為襄陽市西南部位、西部漢江北岸地區、中部漢江東部地區。其中最明顯、面積最大的低溫聚集中心——城市西南方的低溫處為山區，而其他兩處區域則多為森林。可以看到，這些低溫地區的溫度都在30℃以下。由於森林中樹木的遮蔭、吸熱作用，使得森林地表溫度明顯低於其他地區；而對於山區而言，由於氣壓低，空氣稀薄，大氣保溫效果較差；同時高山表面距離地面較遠，無法吸收地面熱量，使得其表面溫度較低。與此同時，二者同樣具有人口密度很低的特點。