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

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

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

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

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

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

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

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

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

Sobrino對大氣校正法、單窗算法和普適性單通道算法進行了比較分析，結(jié)果表明，當(dāng)使用實時大氣廓線數(shù)據(jù)時，大氣校正法的均方根誤差(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。其上游全部分布在祁連山自然保護區(qū)轄區(qū)內(nèi)，在保護區(qū)內(nèi)石羊河上游面積約8 851.5 km2，主要由森林、草原、草甸、荒漠、農(nóng)田、冰川積雪和濕地等自然生態(tài)系統(tǒng)構(gòu)成，是石羊河上游重要水源涵養(yǎng)區(qū)，應(yīng)用遙感技術(shù)進行石羊河上游土地地表溫度反演研究，具有重要科研價值和現(xiàn)實意義。

1.2 影像來源與處理

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

-07-17和2014-07-26，軌道號是 131/034，132/034。Landsat 8衛(wèi)星剛發(fā)射運行不久，TIRS B11波段暫時存在定標(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景影像分別進行輻射定標(biāo)、大氣校正，并進行幾何精校正。校正后利用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ù)進行裁剪。應(yīng)用波段運算得到歸一化植被指數(shù)(NDVI)。因2景影像成像的時間不同，大氣上行輻射、下行輻射以及大氣透過率也不同，應(yīng)分別進行地表溫度的反演，最后將反演的地表溫度影像進行鑲嵌、裁剪，得到研究區(qū)地表溫度反演影像數(shù)據(jù)。

2 研究方法

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

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

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

(1)

式中：ε為地表輻射率；TS為地面真實溫度，單位為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/)輸入影像的成像時間、中心經(jīng)緯度等相關(guān)信息生成，黑體的輻射亮度B(TS)單位為W/(m2·sr· μm) 。

由于NASA官網(wǎng)暫時只能獲取B10波段的參數(shù)，而不能獲取B11波段的參數(shù)，因此統(tǒng)一使用B10波段反演地表溫度。估算出地表真實溫度相同的黑體的輻射亮度B(TS)后，根據(jù)普朗克定律反函數(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軟件直接計算法獲得。植被覆蓋度(Pv)以NDVI數(shù)據(jù)為基礎(chǔ)，采用下式計算。

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

(4)

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

2.2.2 地表比輻射率估算

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

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

2．4 地表溫度反演

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

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

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

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

3 結(jié)果與分析

3.1 地表溫度分級特征

應(yīng)用envi軟件對獲取的研究區(qū)地表溫度進行密度分割(圖1)，研究區(qū)各溫度區(qū)間的統(tǒng)計特征(表1)，研究區(qū)地表溫度在-5.23～50.77℃之間，統(tǒng)計反演結(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個有效取樣點的地表溫度分別與相同點經(jīng)度、緯度、海拔、坡向、坡度作回歸分析，在線性、指數(shù)、對數(shù)、多項式等模型中選擇相關(guān)系數(shù)最大模型作為擬合模型。結(jié)果表明，地表溫度與相同位置的海拔、緯度、坡向相關(guān)關(guān)系極顯著，與相同位置的經(jīng)度、坡度相關(guān)關(guān)系不顯著。并對相關(guān)極顯著的進行模型擬合，擬合結(jié)果見表2，圖2 ～圖4。

表1 研究區(qū)地表溫度密度分割區(qū)間統(tǒng)計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ū)地表溫度—海拔散點圖Fig.2 Scatter diagram of land surface temperature- altitude in the research area

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

圖4 研究區(qū)地表溫度—坡向散點圖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)，地表氣溫低則地表溫度也低，且與實際相吻合。

由圖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)極顯著，表明多項式相關(guān)模型成立。研究區(qū)隨著坡向的變化，山區(qū)陽坡光照時數(shù)長,強度大;陰坡光照時數(shù)短,強度小[10]。因而陽坡的地表溫度、氣溫變幅及水分蒸發(fā)均超過陰坡。陽坡干旱程度超過陰坡，地表植被覆蓋比陰坡差；東、西坡接受的光照時長、強度、土壤水分、植被覆蓋均介于南、北坡之間，造成地表溫度隨坡向出現(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)度、緯度、海拔、坡向、坡度進行相關(guān)、回歸分析。結(jié)果表明，地表溫度隨海拔的升高而降低，呈線性負相關(guān),相關(guān)極顯著；地表溫度隨緯度的增大而升高，呈線性正相關(guān),相關(guān)極顯著；地表溫度隨坡向的變化呈周期性變化趨勢，呈多項式相關(guān),相關(guān)極顯著；地表溫度隨經(jīng)度的增大而降低，相關(guān)不顯著；地表溫度隨坡度的增大而降低，線性相關(guān)不顯著。

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

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

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

[3] 歷華，曾永年，負培東,等.利用多源遙感數(shù)據(jù)反演城市地表溫度[J].遙感學(xué)報，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é)報，2015，19(6)：964-976.

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

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

[9] 毛文婷,王旭紅,祝明英,等.城市地表溫度反演及其與下墊面定量關(guān)系分析—以西安市為例[J].山東農(nóng)業(yè)大學(xué)學(xué)報(自然科學(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è)科技項目(項目編號：2015kj043).

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

10.3969/j.issn.1671-3168.2017.03.005

S716.2；P423

A

1671-3168(2017)03-0017-05