丁國民，汪有奎，邸 華，車宗璽，李碩軍

(甘肅祁連山國家級自然保護(hù)區(qū)管理局，甘肅 張掖 734000)

甘肅石羊河流域上游地表溫度反演及空間特征分析

丁國民，汪有奎，邸 華，車宗璽，李碩軍

(甘肅祁連山國家級自然保護(hù)區(qū)管理局，甘肅 張掖 734000)

采用大氣校正法反演研究區(qū)地表溫度，用地表溫度數(shù)據(jù)分別與經(jīng)度、緯度、海拔、坡向、坡度進(jìn)行相關(guān)、回歸分析，研究地表溫度空間分布特征。以2014年7月的2景Landsat8影像為數(shù)據(jù)源，以甘肅石羊河流域上游為研究對象，應(yīng)用大氣校正法反演了研究區(qū)地表溫度。結(jié)果表明，研究區(qū)地表溫度在-5.23～50.77℃之間，60.94%區(qū)域的地表溫度集中在20～30℃之間；研究區(qū)地表溫度隨海拔的升高而降低，呈線性負(fù)相關(guān),相關(guān)極顯著；地表溫度隨緯度的增大而升高，呈線性正相關(guān),相關(guān)極顯著；地表溫度隨坡向的變化呈周期性變化趨勢，呈多項(xiàng)式相關(guān),相關(guān)極顯著；地表溫度隨經(jīng)度的增大而降低，相關(guān)不顯著；地表溫度隨坡度的增大而降低，線性相關(guān)不顯著。

大氣校正法；地表溫度；反演；空間分布特征；環(huán)境因子；甘肅石羊河流域

陸地表面溫度是地球表面能量平衡和溫室效益的一個重要指標(biāo)。傳統(tǒng)獲取地表溫度的途徑只能依靠人工測量，只能獲取小范圍的點(diǎn)數(shù)據(jù)；遙感數(shù)據(jù)范圍大、全天候等特點(diǎn)為獲取地表溫度數(shù)據(jù)提供了便利手段[1]。地表溫度反演算法主要有3種：大氣校正法(也稱為輻射傳輸方程：Radiative Transfer Equation——RTE)、單通道算法和分裂窗算法。輻射傳輸方程描述了衛(wèi)星的微波輻射計(jì)所觀測到的輻射總強(qiáng)度,不僅有來自地表的輻射,而且還有來自大氣的向上和向下的路徑輻射。這些輻射成分在穿過大氣層到達(dá)遙感器的過程中,還受到大氣層吸收作用的影響而削減。同時(shí),地表和大氣的輻射也在這一過程中產(chǎn)生不可忽略的影響。因此,地表溫度的演算實(shí)際上是一個復(fù)雜的求解問題[2]。

Sobrino對大氣校正法、單窗算法和普適性單通道算法進(jìn)行了比較分析，結(jié)果表明，當(dāng)使用實(shí)時(shí)大氣廓線數(shù)據(jù)時(shí)，大氣校正法的均方根誤差(RMSE)是0.6K，單窗算法和普適性單通道算法的均方根誤差都是0.9K[3，4]；大氣校正法反演地表溫度具有相對較高的精度，故本文基于大氣校正法，利用Landsat8 TIRS反演研究區(qū)地表溫度。

1 研究區(qū)概況與研究數(shù)據(jù)

1.1 研究區(qū)概況

石羊河是發(fā)源于祁連山區(qū)的三大內(nèi)陸河之一，石羊河流域上游的祁連山區(qū)屬高寒半干旱濕潤區(qū)：海拔2 000～5 000 m，年降水量300～600 mm，年蒸發(fā)量700～1 200 mm，干旱指數(shù)1～4。其上游全部分布在祁連山自然保護(hù)區(qū)轄區(qū)內(nèi)，在保護(hù)區(qū)內(nèi)石羊河上游面積約8 851.5 km2，主要由森林、草原、草甸、荒漠、農(nóng)田、冰川積雪和濕地等自然生態(tài)系統(tǒng)構(gòu)成，是石羊河上游重要水源涵養(yǎng)區(qū)，應(yīng)用遙感技術(shù)進(jìn)行石羊河上游土地地表溫度反演研究，具有重要科研價(jià)值和現(xiàn)實(shí)意義。

1.2 影像來源與處理

目前用于反演地表溫度的遙感數(shù)據(jù)主要集中在NOAA、MODIS、TM等數(shù)據(jù)中，這些數(shù)據(jù)空間分辨率都較低，已經(jīng)不能滿足當(dāng)前應(yīng)用的需求[4，5]。本文選取Landsat8遙感數(shù)據(jù)為數(shù)據(jù)源，成像時(shí)間為2014

-07-17和2014-07-26，軌道號是 131/034，132/034。Landsat 8衛(wèi)星剛發(fā)射運(yùn)行不久，TIRS B11波段暫時(shí)存在定標(biāo)不穩(wěn)定性，故本文主要針對 Landsat 8 TIRS B10數(shù)據(jù)討論大氣校正法算法[6]。

1.3 DEM及地形數(shù)據(jù)

采用“地理空間數(shù)據(jù)云”公布的GDEMV2 30M分辨率數(shù)字高程數(shù)據(jù)，通過env 4.8軟件地形分析功能，用DEM數(shù)據(jù)提取研究區(qū)的高程、坡向、坡度數(shù)據(jù)。

1.4 數(shù)據(jù)預(yù)處理

首先將2景影像分別進(jìn)行輻射定標(biāo)、大氣校正，并進(jìn)行幾何精校正。校正后利用envi 4.8軟件自定義坐標(biāo)系，將2景影像、地形數(shù)據(jù)、石羊河上游矢量邊界矢量數(shù)據(jù)的投影和坐標(biāo)系統(tǒng)一為transverse mercator，D-Beijing-1954坐標(biāo)系，應(yīng)用石羊河流域上游矢量數(shù)據(jù)對DEM及地形數(shù)據(jù)進(jìn)行裁剪。應(yīng)用波段運(yùn)算得到歸一化植被指數(shù)(NDVI)。因2景影像成像的時(shí)間不同，大氣上行輻射、下行輻射以及大氣透過率也不同，應(yīng)分別進(jìn)行地表溫度的反演，最后將反演的地表溫度影像進(jìn)行鑲嵌、裁剪，得到研究區(qū)地表溫度反演影像數(shù)據(jù)。

2 研究方法

2.1 LST模型構(gòu)建與參數(shù)計(jì)算

在地—?dú)獾妮椛鋫鬏斨校l(wèi)星接收到的熱紅外輻射能量Lλ包含3部分內(nèi)容：地面真實(shí)輻射經(jīng)大氣衰減之后被衛(wèi)星傳感器接收到的熱輻射能量、大氣的上行輻射亮度L↑、大氣下行輻射亮度L↓(大氣向地面熱輻射)。衛(wèi)星傳感器接收到的熱紅外輻射亮度值Lλ的表達(dá)式(輻射傳輸方程)可寫為：

Lλ= [εB(TS) + (1-ε)L↓]τ + L↑

(1)

式中：ε為地表輻射率；TS為地面真實(shí)溫度，單位為k；τ為大氣在熱紅外波段的透過率。L↓大氣下行輻射亮度和L↑大氣上行輻射亮度，其單位均為W/(m2·sr·μm)。B(TS)為普朗克定律推算得到的黑體在TS的熱輻射亮度。

則溫度為T的黑體在熱紅外波段的輻射亮度B(TS)為：

B(TS) = [Lλ- L↑-τ(1-ε)L↓]/τε

(2)

式中的大氣在熱紅外波段的透過率 τ、大氣上行輻射亮度L↑和大氣下行輻射亮度L↓單位均為W/(m2·sr· μm)，3個參數(shù)可以通過NASA官網(wǎng)(http：//atmcorr．gsfc．nasa．gov/)輸入影像的成像時(shí)間、中心經(jīng)緯度等相關(guān)信息生成，黑體的輻射亮度B(TS)單位為W/(m2·sr· μm) 。

由于NASA官網(wǎng)暫時(shí)只能獲取B10波段的參數(shù)，而不能獲取B11波段的參數(shù)，因此統(tǒng)一使用B10波段反演地表溫度。估算出地表真實(shí)溫度相同的黑體的輻射亮度B(TS)后，根據(jù)普朗克定律反函數(shù)，得出地面真實(shí)溫度，公式為：

TS=K2/ln(K1/B(TS)+ 1)

(3)

式中：TS為傳感器處的地表溫度(K)，B(TS) 為黑體在熱紅外波段的輻射亮度，K1和K2為熱紅外波段的定標(biāo)常數(shù)，對于TIRS 10波段，K1=774.89W/(m2·sr·μm)，K2=1321.08W/(m2·sr·μm)。K1和K2可在下載影像頭文件中查得。

2.2 地表比輻射率及參數(shù)估算2.2.1 NDVI及植被覆蓋度估算

NDVI指數(shù)采用envi軟件直接計(jì)算法獲得。植被覆蓋度(Pv)以NDVI數(shù)據(jù)為基礎(chǔ)，采用下式計(jì)算。

Pv=[(NDVI- NDVISoil)/(NDVIVeg-NDVISoil)]

(4)

式中：NDVI為歸一化植被指數(shù)，NDVISoil為完全是裸土或無植被覆蓋區(qū)域的NDVI值，NDVIVeg則代表完全被植被所覆蓋的像元的NDVI值，即純植被像元的NDVI值，取經(jīng)驗(yàn)值NDVIVeg= 0.70和NDVISoil= 0.05，即當(dāng)某個像元的NDVI大于0.70時(shí)，Pv取值為1；當(dāng)NDVI小于0.05，Pv取值為0。

2.2.2 地表比輻射率估算

地表比輻射率是地表溫度反演中不可或缺的重要參數(shù)，目前對地表比輻射率的估計(jì)方法主要是基于NDVI的地表比輻射率法。本文采用覃志豪等提出的地表比輻射率計(jì)算方法，先將地表分成水體、自然表面和自然、人工混合表面，分別針對3種地表類型計(jì)算地表比輻射率[9]：

ε=0.995 (NDVI≤NDVIv) (水體、冰雪覆蓋區(qū))

ε= 0.9589 + 0.086Pv-0.0671Pv2(NDVIv

ε=0.9625 + 0.0614Pv-0.0461Pv2(NDVI>NDVIs)(自然表面)

式中：Pv是植被覆蓋度；NDVIs 為完全被裸土或無植被覆蓋區(qū)域的NDVI 值，NDVIv 則代表完全被植被所覆蓋的像元的NDVI 值，即純植被像元的NDVI值。

2.3 黑體輻射亮度計(jì)算

大氣校正法進(jìn)行地表溫度的反演，其中涉及到大氣上行輻射、下行輻射以及大氣透過率數(shù)據(jù)，這些數(shù)據(jù)可以在NASA官網(wǎng)(http://atmcorr.gsfc.nasa.gov/)中輸入成像的時(shí)間及其中心經(jīng)度，則會得到所需要的參數(shù)。將相應(yīng)的參數(shù)帶入公式2，即可得到同溫度下的黑體輻射亮度圖像。

2．4 地表溫度反演

將Landsat 8同溫度下的TIRS 10黑體輻射亮度圖像代入公式(3)，即可獲得研究區(qū)地表溫度。應(yīng)用envi軟件對地表溫度進(jìn)行分級,得到地表溫度分級圖(圖1)。

圖1 研究區(qū)地表溫度密度分級Fig.1 Classification of land surface temperature and density in the research area

2.5 地表溫度空間分布特征分析方法

應(yīng)用ArcGIS10.2軟件，對鑲嵌、裁剪后的地表溫度反演圖進(jìn)行魚網(wǎng)取樣，取樣的間隔為0.8 km，獲得魚網(wǎng)點(diǎn)經(jīng)度、緯度、地表溫度數(shù)據(jù)。用同樣的方法應(yīng)用ArcGIS10.2對海拔、坡度和坡向影像圖進(jìn)行魚網(wǎng)采樣，取樣精度小于1個像素，獲得相同經(jīng)度、緯度下的海拔、坡向、坡度數(shù)據(jù)，用于地表溫度空間分布特征分析。用地表溫度數(shù)據(jù)分別與經(jīng)度、緯度、海拔、坡向、坡度進(jìn)行相關(guān)、回歸分析，得到地表溫度與經(jīng)度、緯度、海拔、坡向、坡度的回歸模型。

3 結(jié)果與分析

3.1 地表溫度分級特征

應(yīng)用envi軟件對獲取的研究區(qū)地表溫度進(jìn)行密度分割(圖1)，研究區(qū)各溫度區(qū)間的統(tǒng)計(jì)特征(表1)，研究區(qū)地表溫度在-5.23～50.77℃之間，統(tǒng)計(jì)反演結(jié)果得出0.84%區(qū)域的地表溫度在10℃以下，1.98%區(qū)域的地表溫度在10～20℃之間，60.94%區(qū)域的地表溫度集中在20～30℃之間，29.95%區(qū)域的地表溫度在30～40℃之間，6.28%區(qū)域的地表溫度在40℃以上。

3.2 地表溫度與環(huán)境因子的相關(guān)模型

采用133個有效取樣點(diǎn)的地表溫度分別與相同點(diǎn)經(jīng)度、緯度、海拔、坡向、坡度作回歸分析，在線性、指數(shù)、對數(shù)、多項(xiàng)式等模型中選擇相關(guān)系數(shù)最大模型作為擬合模型。結(jié)果表明，地表溫度與相同位置的海拔、緯度、坡向相關(guān)關(guān)系極顯著，與相同位置的經(jīng)度、坡度相關(guān)關(guān)系不顯著。并對相關(guān)極顯著的進(jìn)行模型擬合，擬合結(jié)果見表2，圖2 ～圖4。

表1 研究區(qū)地表溫度密度分割區(qū)間統(tǒng)計(jì)Tab.1 Interval partition of land surface temperature and density in the research area

表2 地表溫度與生態(tài)因子的相關(guān)模型Tab.2 Model of land surface temperature and ecological factor

圖2 研究區(qū)地表溫度—海拔散點(diǎn)圖Fig.2 Scatter diagram of land surface temperature- altitude in the research area

圖3 研究區(qū)地表溫度—緯度散點(diǎn)圖Fig.3 Scatter diagram of land surface temperature- latitude in the research area

圖4 研究區(qū)地表溫度—坡向散點(diǎn)圖Fig4 Scatterdiagramoflandsurfacetemperature-aspectintheresearcharea

地表溫度隨經(jīng)度的增加而降低，相關(guān)系數(shù)r=0.151 0

由圖2可知，研究區(qū)地表溫度隨海拔的升高而降低，相關(guān)系數(shù)r=0.3191>r0.01=0.222,相關(guān)極顯著，表明線性相關(guān)模型成立，地表溫度的高低與地表氣溫相關(guān)，地表氣溫低則地表溫度也低，且與實(shí)際相吻合。

由圖3可知，研究區(qū)地表溫度隨緯度的增大而升高，相關(guān)系數(shù)r=0.413 9>r0.01=0.222,相關(guān)極顯著，表明線性相關(guān)模型成立。研究區(qū)隨著緯度增加、海拔降低，地表溫度逐步升高，這是海拔、緯度變化雙重作用所致。祁連山緯度對氣溫的影響比較特殊，一般緯度增加，氣溫會降低，其中由于海拔高度變化的影響，祁連山區(qū)出現(xiàn)了緯度增加，氣溫升高的“逆反”現(xiàn)象。

由圖4可知，研究區(qū)地表溫度隨坡向的變化呈現(xiàn)周期性變化趨勢，相關(guān)系數(shù)r=0.319 7>r0.01=0.222,相關(guān)極顯著，表明多項(xiàng)式相關(guān)模型成立。研究區(qū)隨著坡向的變化，山區(qū)陽坡光照時(shí)數(shù)長,強(qiáng)度大;陰坡光照時(shí)數(shù)短,強(qiáng)度小[10]。因而陽坡的地表溫度、氣溫變幅及水分蒸發(fā)均超過陰坡。陽坡干旱程度超過陰坡，地表植被覆蓋比陰坡差；東、西坡接受的光照時(shí)長、強(qiáng)度、土壤水分、植被覆蓋均介于南、北坡之間，造成地表溫度隨坡向出現(xiàn)周期性的變化。

4 結(jié)語

利用研究區(qū)域Landsat 8數(shù)據(jù)(使用B10波段)結(jié)合輻射傳導(dǎo)方程法得到石羊河流域上游地表溫度，反演結(jié)果表明，研究區(qū)60.94%的區(qū)域處于中溫度區(qū)域。高溫區(qū)、低溫區(qū)所占的比例分別為6.28%和0.84%，所占的面積均較低。

用地表溫度數(shù)據(jù)分別與經(jīng)度、緯度、海拔、坡向、坡度進(jìn)行相關(guān)、回歸分析。結(jié)果表明，地表溫度隨海拔的升高而降低，呈線性負(fù)相關(guān),相關(guān)極顯著；地表溫度隨緯度的增大而升高，呈線性正相關(guān),相關(guān)極顯著；地表溫度隨坡向的變化呈周期性變化趨勢，呈多項(xiàng)式相關(guān),相關(guān)極顯著；地表溫度隨經(jīng)度的增大而降低，相關(guān)不顯著；地表溫度隨坡度的增大而降低，線性相關(guān)不顯著。

由于沒有同步實(shí)測數(shù)據(jù)，無法對反演的地表溫度作出定量的評價(jià)，但研究得出地表溫度與地形因子的經(jīng)度、緯度、海拔、坡向、坡度之間的相關(guān)關(guān)系和回歸模型及地表溫度空間分布特征，仍有一定的參考價(jià)值。

[1] 蔣大林，匡鴻海，曹小峰，等.基于Landsat 8的地表溫度反演算法研究—以滇池流域?yàn)槔齕J].遙感技術(shù)與應(yīng)用，2015,30(3)：448-454.

[2] 宋金紅,陳圣波,包書新,等.地表溫度反演的數(shù)值方法[J].吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2007(增刊)：198-201.

[3] 歷華，曾永年，負(fù)培東,等.利用多源遙感數(shù)據(jù)反演城市地表溫度[J].遙感學(xué)報(bào)，2007，11(6)：891-898.

[4] 樊輝．基于Landsat TM的城市熱島效應(yīng)與地表特征參數(shù)穩(wěn)健關(guān)系模型[J]．國土資源遙感，2008，77(3)：52-57.

[5] 鄒蒲，李婷，鄒娟.基于Landsat8衛(wèi)星數(shù)據(jù)的城市熱島效應(yīng)分析—以湘潭市為例[J]．環(huán)球人文地理(評論版)，2015(3)：32-34.

[6] 胡德勇，喬琨，王興玲,等.單窗算法結(jié)合 Landsat 8 熱紅外數(shù)據(jù)反演地表溫度[J].遙感學(xué)報(bào)，2015，19(6)：964-976.

[7] 吳志剛，江滔，樊艷磊,等.基于Landsat 8數(shù)據(jù)的地表溫度反演及分析研究—以武漢市為例[J].工程地球物理學(xué)報(bào)，2016，13(1)：135-142.

[8] 游絢，晏路明.基于ETM+ 影像輻射傳導(dǎo)方程算法的地表溫度反演[J].科技情報(bào)開發(fā)與經(jīng)濟(jì)，2009，19(27)：134-136．

[9] 毛文婷,王旭紅,祝明英,等.城市地表溫度反演及其與下墊面定量關(guān)系分析—以西安市為例[J].山東農(nóng)業(yè)大學(xué)學(xué)報(bào)(自然科學(xué)版),2015,46(5):708-714.

[10] 徐化成.景觀生態(tài)學(xué)[M].北京：中國林業(yè)出版社,1995：1-83.

Inversion and Spatial Characteristics of Land Surface Temperatures inthe Upper Shiyang River Basin of Gansu

DING Guomin, WANG Youkui, DI Hua, CHE Zongxi LI Shuojun

(Gansu Qilian Mountains National Natural Reserve, Zhangye, Gansu 734000,China)

In order to study the spatial distribution characteristics of land surface temperatures, we usually study the land surface temperature by method of atmospheric correction, and make correlation and regression analysis based on the land surface temperature data with longitude, latitude, altitude, aspect, and gradient respectively. Taking 2 Landsat 8 images of July 2014 as the data resource, the upper Shiyang river basin of Gansu as the research object, surface temperatures were displayed by the method of atmospheric correction. The results showed that land surface temperatures in the research area ranged from -5.23 to 50.77, among which land surface temperatures of 60.94% areas were between 20 and 30 degrees; the land surface temperature decreased with the increase of altitude, presenting the remarkable linear negative correlation, heightened with the increase of latitude, presenting the remarkable linear positive correlation, showed periodic variation trend with the change of aspect, markedly presenting the polynomial correlation, decreased with the increase of longitude, inapparently presenting the correlation; and decreased with the increase of gradient, indistinctively presenting the linear correlation.

atmospheric correction; land surface temperatures; inversion; spatial distribution characteristics; environmental factors ;Shiyang river basin of Gansu

2017-03-30.

2015甘肅省林業(yè)科技項(xiàng)目(項(xiàng)目編號：2015kj043).

丁國民(1966-)，男，甘肅人，高級工程師.從事保護(hù)區(qū)生態(tài)保護(hù)及科研工作.

10.3969/j.issn.1671-3168.2017.03.005

S716.2；P423

A

1671-3168(2017)03-0017-05