大興安嶺森林地表溫度的遙感估算及分析研究

孫舒婷，范文義，張智超

(東北林業大學林學院，哈爾濱150040)

地表能量交換信息的獲取是監測區域森林資源環境變化的一個重要環節。獲取區域地表溫度的空間差異并進行分析，評價森林對溫度的調節作用是對森林生態系統生態服務功能研究的重要內容［1－3］。目前，對森林調節溫度的生態服務功能的研究多數在點尺度上利用氣象觀測點的數據進行研究，對區域尺度森林的溫度調節作用定量評價研究需要尋求新的手段。

熱紅外遙感技術的飛速發展為快速地獲取區域地表溫度差異信息提供了新途徑。美國NASA的地球觀測計劃提出了陸面溫度反演精度優于1K，海面溫度反演精于優于0.3K的目標，認為這個精度能夠對遙感應用產生實質的推動。Landsat陸地衛星的TM遙感影像數據具有地面分辨率高的特點，TM6熱紅外波段在分析地表溫度區域差異已經有了廣泛的應用。對于單波段熱紅外較成熟的算法大致分為三類，輻射傳輸方程法、單波段算法和普適性單窗算法，其產品廣泛應用于城市熱島效應等相關研究［4－6］，但對于森林生態服務功能的調節溫度效應方面研究還較少。

本文利用TM的熱紅外波段、植被指數及大興安嶺二類調查數據和氣象數據等，對大興安嶺地區森林對溫度調節的生態服務功能進行量化研究，為科學地評估森林生態系統調節溫度的生態服務功能提供新的技術手段和科學的數據支持，初步探究適用于森林地區的地表溫度反演方法。

1 研究區概況及數據收集處理

大興安嶺位于黑龍江省、內蒙古自治區東北部(東經 121°12'～ 127°00';北緯 50°10'～ 53°33')，是中國最重要的林業基地之一。它北起黑龍江畔，南至西林木河上游谷地，全長1 200 km，寬200～300 km，海拔1 100～1 400 m，是中國面積最大的林區，木材貯量占中國的一半。大興安嶺的林地有730萬hm2，森林覆蓋率約74%，林層結構復雜，林分優勢樹種以興安落葉松、樟子松、白樺為主。

本研究所采用大興安嶺地區Landsat TM影像、氣象站數據和二類調查數據。Landsat－5軌道高度705km，軌道傾角98.22°，掃描帶寬為185km，每個波段具體參數見表1。

表1 Landsat－5主要參數Tab.1 Major parameters of Landsat－5

由于每景TM數據的獲取時間不一致，特別是不同軌道之間的數據獲取時間差別更大，用不同景TM數據估算溫度后沒法比較，所以本文選取其中比較有代表性的一景TM影像來進行說明分析。該影像的獲取時間為2007年7月23號，中心坐標為124.99E，52.796N。氣象站數據在中國氣象科學數據共享服務網下載。二類調查數據是2006年大興安嶺森林資源二類調查數據，包括立地類型和林相圖等數據。

對所選取的TM影像進行輻射校正和幾何精校正，結合地形圖進行幾何精校正處理，得到結果精度控制在0.5個像元內。按公式(1)求得植被指數，結果如圖1所示。

圖1 NDVIFig.1 NDVI

式中:TM4和TM3是Landsat TM影像的第4波段和第3波段的DN值。

因為只需要運用NDVI這個歸一化的指數，所以不需要再進行大氣校正［7］。用真彩色合成可以預覽整幅區域的概況，如圖2所示，除森林覆蓋區域之外，該區域內有一些建筑用地、裸地和水面。

圖2 真彩色合成(RGB=bands3，2，1)Fig.2 Real color combination(RGB=bands3，2，1)

2 研究方法

2.1 TM6地表溫度反演方法

對于植被茂密的地表，遙感反演所得到的地表溫度是指植被葉冠層的表面溫度。對于稀疏的地表，地表溫度是地面、植被葉冠等溫度的混合平均值。因此，地表的非同質性使地表溫度的遙感反演成為一個復雜的問題，目前遙感能獲取的仍然是像元的平均溫度，所以TM數據的高空間分辨率就體現了它的優勢。

2.1.1 輻射傳輸方程法

如果能夠獲取大氣的溫、濕度廓線，使用大氣輻射傳輸方程［8］，可以模擬出公式(2)中的3個大氣參數:上行輻射亮度、下行輻射亮度和透射率，假設地表發射率已知，可由普朗克函數求逆得到地表溫度。

式中:Li為衛星高度上傳感器測得的輻射強度(W/m2·sr·μm4);εi為地表輻射率;Bi(Ts)為由Plank定律推導得到的黑體熱輻射強度，其中Ts為地表溫度 (K);Lu，i和 Ld，i分別為大氣上行輻射亮度和大氣下行輻射亮度;τi為大氣透射率。

這種算法實際操作起來非常困難，除計算過程復雜外，大氣模擬所需要的數據實時大氣氣溫、濕度廓線數據很難獲取。

2.1.2 單波段算法

Qin和Karnieli［9］提出了針對TM的單窗算法。根據熱輻射傳輸方程使用中值定理［10］，引入大氣平均作用溫度Ta來近似表達大氣上行輻射亮度和下行輻射亮度。假設大氣向上的平均作用溫度和向下的平均作用溫度相等，并在常溫下對普朗克函數線性近似，得到地表溫度的表達式:

式中:Ts是地表溫度;T6是TM第六波段的亮度溫度;Ta是大氣平均作用溫度;a=－67.355 351;b=0.458 606;C和D為中間變量:C=ε6τ6;D=(1－τ6) ［1+τ6(1－ ε6)］;ε6和 τ6分別是第六波段的地表發射率和大氣透射率。

如果TM6圖象的溫度變化范圍較窄［4］，還可提高估計誤差。例如，對于0～30℃，取 a6=－60.326 3和 b6=0.434 36。

該算法僅需要知道3個參數:地表發射率，大氣透射率和平均作用溫度。大氣透過率和平均作用溫度可由大氣溫、濕度輪廓線或氣象站點的觀測數據估算。

2.1.3 普適性單通道算法

Jiménez－Mu?oa和 Sobrino提出了一個普適性的單通道算法［11－12］，該算法可以針對任何一種熱紅外數據反演地表溫度，同樣適用于TM6數據。算法中地表溫度表示為:

與單波段算法相比，該算法更為簡單，所需的輸入參數除地表發射率外僅需要大氣水含量。

對比以上3種算法，輻射傳輸方程法需要實時大氣氣溫、濕度廓線數據，這些數據不易獲取且誤差相對較大;單波段算法運用線性近似的方法，在地表溫度差異較大時也會增大誤差;普適性單窗算法需要的數據易獲取操作簡單方便并且誤差較低。這與Sobrino等人在2004年的研究結果吻合［7］，當使用實時大氣廓線數據時，輻射傳輸方程的均方根誤差(RMSE)為0.6K，單窗算法和普適性單通道算法的RMSE均為0.9K;當沒有實時大氣廓線數據時，輻射傳輸方程法不再適用，而單窗算法和普適性單通道算法的RESE分別為2K和0.9K。基于以上3種方法和分析本研究將采用普適性單通道算法。

2.2 地表輻射參數估算

對于自然界絕大多數地表，發射率具有方向性，在地表發射率的遙感反演中，通常假設地表為朗伯體，忽略其方向性。考慮到大興安嶺地區的地貌，從衛星像元尺度看，大體可分為3種類型:水面、建筑用地和植被覆蓋的自然表面。

水體在熱波段范圍內的比輻射率很高，接近于黑體［13］，可以用εw=0.995來進行估計。對于建筑用地表面而言，在TM6波段范圍內，其比輻射率一般在0.960～0.980之間變動，這里用εm=0.970代替。實際上，組成自然表面的像元可以簡單地看作是由不同比例的植被葉冠和裸土組成。一般通過下式來估算自然表面的混合地表比輻射率［14］。

式中:Pv是植被占混合像元的比例;Rv和Rs分別是植被和裸土的溫度比率;εv和εs分別是植被和裸土的地表比輻射率。

根據Labed and Stoll和Humes等人的測量結果［15－16］，植被的比輻射率可以用 εv=0.986 進行初步估計［17－18］;裸土在TM6波段區的比輻射率用棕壤土、粘質土、沙質土和沙壤土的平均值εs=0.972 15代替。

植被占混合像元構成比率Pv的估算大致可以分為3種:經驗模型法、植被指數法和像元分解模型法。經驗模型法的精度高但對實測數據有依賴性適用于小區域精確研究;植被指數法的精度較低，對實測數據依賴性較小，適用于大范圍粗略估計;像元分解模型法精度隨著分辨率而不同，不依賴實測數據，目前應用范圍最為廣泛。結合本研究的區域和數據，采用像元分解模型法中簡單實用的像元二分模型對植被構成比率Pv做估測。假設像元只由植被與非植被兩種地表覆蓋類型，它們各自的面積在像元中所占的像元百分比即為Pv，表達式為:

式中:NDVI是歸一化植被指數;NDVIs和NDVIv分別是裸土和植被的NDVI值。

NDVI值越大，越接近于完全植被覆蓋;NDVI值越小，越接近于完全裸土覆蓋;當NDVI介于植被和裸土之間時，表明有一定比例的植被覆蓋和一定比例的裸土。由于NDVI是經過歸一化的植被指數，大氣影響對于NDVI影響不大，NDVI誤差引起的地表溫度誤差極小(﹤0.1K)［7］，因此，不用進行大氣校正，可以直接用TM3和TM4的DN值來計算。

3 研究結果

3.1 地表溫度的反演

先求得星上亮度溫度和NDVI如圖1和圖3所示，大氣水含量ω由氣象站數據獲得;根據大興安嶺地區NDVI分布圖和二類調查數據的立地類型求得，明顯的濃密植被區的平均NDVI值作為NDVIv值，即NDVIv=0.71;明顯的裸土區的平均NDVI值作為 NDVIs值;即 NDVIs=0.05。當 NDVI＞NDVIV時，PV=1;當 NDVI＜NDVIS時，PV=0;NDVIV≤NDVI≤NDVIS時，該像元認為是由裸土的植被組成的混合像元，此時用公式(10)～(13)計算。運用更精確的普適性單窗算法，得到該幅影像的溫度分布，如圖4所示，將溫度分布圖分為4個梯度:20℃以下;20～25℃;25～30℃;30℃以上，得到溫度梯度圖，如圖5所示。

圖3 星上亮度溫度Fig.3 Satellite radiance temperature

圖4 溫度分布Fig.4 Temperature distribution

圖5 溫度梯度圖Fig.5 Temperature gradients

由于無法直接獲取地表溫度來驗證精度，本文用氣象站獲取的氣溫數據進行檢驗分析。由于2007年7月23日大興安嶺地區氣象數據中最高溫度和最低溫度的觀測站點較少，通過大區域插值并裁剪獲得與TM影像范圍相同的當日最高氣溫和最低氣溫，與經過幾何精校正的地表溫度產品進行分析。將反演得到的地表溫度柵格值按溫度從低到高排列，加入對應的最高溫度與最低溫度，統計結果如圖6所示，從圖中可以看到，反演出的地表溫度介于當日最高溫度和最低溫度之間，Landsat陸地衛星過境的時間在10點40分左右，該時間的地表溫度介于最高溫度和最低溫度之間，反演出的地表溫度曲線平滑，接近于真實溫度值，提供了比氣象站更理想的空間異質度信息。

圖6 溫度趨勢圖Fig.6 The tendency of temperature variation

3.2 森林調節溫度的生態服務功能評價分析

因為二類調查數據不能覆蓋整景Landsat TM影像，所以運用監督分類法與二類調查的立地類型數據相結合將該影像分為四種地表類型:水面、林地、建筑用地和裸地，如圖7所示，其中藍色為水面、綠色為林地、黃色為建筑用地、紅色色為裸地。水面的占地面積887 364.0 m2，占整幅圖像的2.8%，平均溫度21.5℃，最低溫度16.6℃，最高溫度31.4℃;林地的占地面積為24 799 053.6 m2，占整幅圖像的77.3%，平均溫度為24.8℃，最低溫度18.1℃，最高溫度35.4℃;建筑用地的占地面積2 442 925.8 m2，占整幅圖像的7.6%，平均溫度32.7℃，最低溫度25.0℃，最高溫度44.1℃;裸地的占地面積3 931 419.6 m2，占整幅圖像的12.2%，平均溫度29.3℃，最低溫度18.2℃，最高溫度42.5℃。可以清楚地看到森林在夏季有降溫作用，森林平均溫度比城鎮的建筑用地降低7.9℃。

圖7 立地類型分布圖Fig.7 Distribution of different stand types

4 結論與討論

本研究利用遙感數據、氣象站數據和大興安嶺森林資源二類調查數據，得到科學的地表溫度產品，對定量評價森林生態服務功能提供了科學的依據和新的技術手段，初步探索適用于大興安嶺地區的地表溫度反演方法，得到地表溫度產品，得出森林在夏季有降溫作用，森林平均溫度比城鎮的建筑用地平均溫度降低7.9℃的結論，但對于森林調節溫度的原理和機制沒有涉及，對于大范圍遙感數據的時間不統一問題沒有解決。

通常所研究的森林調節溫度的生態服務中的溫度是指氣溫，遙感反演的溫度是地表溫度，雖然氣溫與地溫存在極強的相關性，但也有所差別，目前沒有成熟的可以借鑒的轉換公式。大興安嶺地區海拔分布在1 100～1 400 m，屬淺山丘陵地帶，本研究沒有考慮高程對地表溫度的影像。TM5溫度產品的分辨率在120m，城鎮邊緣存在混合像元，森林調節溫度的能力在一定程度上被削弱。森林調節溫度的生態服務功能不僅體現在夏季的降溫作用，在冬季同樣具有保溫作用由于數據的限制沒能進行分析。森林的調節溫度作用程度的高低與樹種、郁閉度、地形、時間和環境等很多要素都有關系，但本研究只選取一個時間點定量的進行分析，詳細的數據分析仍然需要實際的地面觀測數據。

［1］李文華，李 芬，李世東，等.森林生態效益補償的研究現狀與展望［J］.自然資源學報，2006，21(5):677 －688.

［2］彭 建，王仰麟，陳燕飛，等.城市生態系統服務功能價值評估初探－以深圳市例［J］.北京大學學報:自然科學版，2005，41(4):594－604.

［3］ Brack C L.Pollution mitigation and carbon sequestration by an urban forest［J］.Envir Poll，2002，116:196 － 200.

［4］ 覃志豪，Zhang Minghua，Karnieli A，等.用陸地衛星 TM6數據演算地表溫度的單窗算法［J］.地理學報，2001，56(4):456－466.

［5］毛克彪，覃志豪，施建成.用MODIS影像和劈窗算法反演山東半島的地表溫度［J］.中國礦業大學學報:自然科學版，2005，34(1):46－50.

［6］張順謙，周長艷.成都市晴天熱島效應的時空分布特征與成因［J］.應用生態學報，2013，24(7):1962 －1968.

［7］ Sobrino J A，Jimnez－Munoza J C，Paolini L.Land surface temperature retrieval from Landsat TM 5［J］.Remote Sens.Environ.，2004，90(4):434－440.

［8］陳良富，莊家禮，徐希孺.熱紅外遙感中通道信息相關性及其對陸面溫度反演的影響［J］.科學通報，1999，44(19):2122 －2127.

［9］ Qin Z，Karnieli A.Mono－window algorithm for retrieving land surface temperature from Landsat TM data and its application to the Israel－Egype border region［J］.International Journal of Remote Sensing，2001，22:3719 －3746.

［10］ McMillin L M.Estimation of sea surface temperature from two infrared window measurements with different absorption［J］.Journal of Geophysical Research，1975，20:11587 －11601.

［11］ Jimenez－Munoz J C，Sobrino J A.A generalized single－channel method for retrieving land surface temperature from remote sensing data［J］.Journal of Geophysical Research，2003，108(D22)，4688.doi:10.1029/2003JD003480.

［12］ Jimenez－Munoz J C，Cristobal J，et al.Revision of the single－channel algorithm for land surface temperature retrieval from Landsat Thermal－Infrared data［J］.IEEE Transactions on Geoscience and Remote Sensing，2009，47:399 －349.

［13］ Cristòbal J，Jimenez－Munoz J C，Sobrino J A，et al.Improvements in land surface temperature retrieval from the Landsat series thermal band using water vapor and air temperature［J］.Journal of Geophysical Research， 2009， 114 (D8). doi:10.1029/2008JD010616.

［14］覃志豪，李文娟，徐 斌，等.利用 Landsat TM6反演地表溫度所需地表輻射率參數的估計方法［J］.海洋科學進展，2004，2(zl):129－137.

［15］ Labed J，Stoll M P.Spatial variability of land surface emissivity in the thermal infrared band:spectral signature and effective surface temperature［J］.Remote Sens.Environ.，1991，38:1 －17.

［16］ Humes K S，Kustas W P，Moran M S，et al.Variability of emissivity and surface temperature over a sparsely vegetated surface［J］.Water Resources Res.，1994，30:1299 －1310.

［17］ Sobrino J A，Raissouni N，Li Z L.A comparative study of land surface emissivity retrieval from NOAA data［J］.Remote Sens.Environ.，2001，75:256 －266.

［18］張淑芬，邢艷秋，艾合買提江·阿不都艾尼，等.基于TM遙感影像的森林類型分類方法比較［J］.森林工程，2014，30(1):18－20.