孫舒婷，范文義，張智超

(東北林業(yè)大學(xué)林學(xué)院，哈爾濱150040)

地表能量交換信息的獲取是監(jiān)測區(qū)域森林資源環(huán)境變化的一個重要環(huán)節(jié)。獲取區(qū)域地表溫度的空間差異并進(jìn)行分析，評價森林對溫度的調(diào)節(jié)作用是對森林生態(tài)系統(tǒng)生態(tài)服務(wù)功能研究的重要內(nèi)容［1－3］。目前，對森林調(diào)節(jié)溫度的生態(tài)服務(wù)功能的研究多數(shù)在點(diǎn)尺度上利用氣象觀測點(diǎn)的數(shù)據(jù)進(jìn)行研究，對區(qū)域尺度森林的溫度調(diào)節(jié)作用定量評價研究需要尋求新的手段。

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

本文利用TM的熱紅外波段、植被指數(shù)及大興安嶺二類調(diào)查數(shù)據(jù)和氣象數(shù)據(jù)等，對大興安嶺地區(qū)森林對溫度調(diào)節(jié)的生態(tài)服務(wù)功能進(jìn)行量化研究，為科學(xué)地評估森林生態(tài)系統(tǒng)調(diào)節(jié)溫度的生態(tài)服務(wù)功能提供新的技術(shù)手段和科學(xué)的數(shù)據(jù)支持，初步探究適用于森林地區(qū)的地表溫度反演方法。

1 研究區(qū)概況及數(shù)據(jù)收集處理

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

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

表1 Landsat－5主要參數(shù)Tab.1 Major parameters of Landsat－5

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

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

圖1 NDVIFig.1 NDVI

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

因?yàn)橹恍枰\(yùn)用NDVI這個歸一化的指數(shù)，所以不需要再進(jìn)行大氣校正［7］。用真彩色合成可以預(yù)覽整幅區(qū)域的概況，如圖2所示，除森林覆蓋區(qū)域之外，該區(qū)域內(nèi)有一些建筑用地、裸地和水面。

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

2 研究方法

2.1 TM6地表溫度反演方法

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

2.1.1 輻射傳輸方程法

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

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

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

2.1.2 單波段算法

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

式中: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分別是第六波段的地表發(fā)射率和大氣透射率。

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

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

2.1.3 普適性單通道算法

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

與單波段算法相比，該算法更為簡單，所需的輸入?yún)?shù)除地表發(fā)射率外僅需要大氣水含量。

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

2.2 地表輻射參數(shù)估算

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

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

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

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

植被占混合像元構(gòu)成比率Pv的估算大致可以分為3種:經(jīng)驗(yàn)?zāi)Ｐ头ā⒅脖恢笖?shù)法和像元分解模型法。經(jīng)驗(yàn)?zāi)Ｐ头ǖ木雀叩珜?shí)測數(shù)據(jù)有依賴性適用于小區(qū)域精確研究;植被指數(shù)法的精度較低，對實(shí)測數(shù)據(jù)依賴性較小，適用于大范圍粗略估計(jì);像元分解模型法精度隨著分辨率而不同，不依賴實(shí)測數(shù)據(jù)，目前應(yīng)用范圍最為廣泛。結(jié)合本研究的區(qū)域和數(shù)據(jù)，采用像元分解模型法中簡單實(shí)用的像元二分模型對植被構(gòu)成比率Pv做估測。假設(shè)像元只由植被與非植被兩種地表覆蓋類型，它們各自的面積在像元中所占的像元百分比即為Pv，表達(dá)式為:

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

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

3 研究結(jié)果

3.1 地表溫度的反演

先求得星上亮度溫度和NDVI如圖1和圖3所示，大氣水含量ω由氣象站數(shù)據(jù)獲得;根據(jù)大興安嶺地區(qū)NDVI分布圖和二類調(diào)查數(shù)據(jù)的立地類型求得，明顯的濃密植被區(qū)的平均NDVI值作為NDVIv值，即NDVIv=0.71;明顯的裸土區(qū)的平均NDVI值作為 NDVIs值;即 NDVIs=0.05。當(dāng) NDVI＞NDVIV時，PV=1;當(dāng) NDVI＜NDVIS時，PV=0;NDVIV≤NDVI≤NDVIS時，該像元認(rèn)為是由裸土的植被組成的混合像元，此時用公式(10)～(13)計(jì)算。運(yùn)用更精確的普適性單窗算法，得到該幅影像的溫度分布，如圖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

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

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

3.2 森林調(diào)節(jié)溫度的生態(tài)服務(wù)功能評價分析

因?yàn)槎愓{(diào)查數(shù)據(jù)不能覆蓋整景Landsat TM影像，所以運(yùn)用監(jiān)督分類法與二類調(diào)查的立地類型數(shù)據(jù)相結(jié)合將該影像分為四種地表類型:水面、林地、建筑用地和裸地，如圖7所示，其中藍(lán)色為水面、綠色為林地、黃色為建筑用地、紅色色為裸地。水面的占地面積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℃?？梢郧宄乜吹缴衷谙募居薪禍刈饔?，森林平均溫度比城鎮(zhèn)的建筑用地降低7.9℃。

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

4 結(jié)論與討論

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

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

［1］李文華，李 芬，李世東，等.森林生態(tài)效益補(bǔ)償?shù)难芯楷F(xiàn)狀與展望［J］.自然資源學(xué)報，2006，21(5):677 －688.

［2］彭 建，王仰麟，陳燕飛，等.城市生態(tài)系統(tǒng)服務(wù)功能價值評估初探－以深圳市例［J］.北京大學(xué)學(xué)報:自然科學(xué)版，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，等.用陸地衛(wèi)星 TM6數(shù)據(jù)演算地表溫度的單窗算法［J］.地理學(xué)報，2001，56(4):456－466.

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

［6］張順謙，周長艷.成都市晴天熱島效應(yīng)的時空分布特征與成因［J］.應(yīng)用生態(tài)學(xué)報，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］陳良富，莊家禮，徐希孺.熱紅外遙感中通道信息相關(guān)性及其對陸面溫度反演的影響［J］.科學(xué)通報，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反演地表溫度所需地表輻射率參數(shù)的估計(jì)方法［J］.海洋科學(xué)進(jìn)展，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.