唐太斌， 周 保， 金曉媚， 魏賽拉加， 馬 濤， 張永艷

（1.青海省地質(zhì)環(huán)境監(jiān)測總站，青海 西寧 810008；2.中國地質(zhì)大學(xué)（北京）水資源與環(huán)境學(xué)院，北京 100083）

黃河源區(qū)作為黃河的發(fā)源地，是中華水塔的重要組成部分。源區(qū)生態(tài)環(huán)境敏感、脆弱，近幾十年來，源區(qū)多年凍土已發(fā)生大面積退化并日趨嚴(yán)重［1］，植被范圍逐漸減少，草場退化、水土流失嚴(yán)重等生態(tài)環(huán)境問題日趨嚴(yán)重。地表溫度是地表、大氣及太陽之間的熱量交換，不少學(xué)者將其視為監(jiān)測地球環(huán)境動態(tài)變化的重要指標(biāo)之一。有學(xué)者將地表溫度應(yīng)用于災(zāi)害監(jiān)測，如利用熱紅外異常來研究地表溫度與地震發(fā)生的關(guān)系［1］。研究黃河源區(qū)地表溫度是深入研究三江源地區(qū)能量循環(huán)、環(huán)境保護(hù)的基礎(chǔ)，可為黃河源區(qū)的生態(tài)保護(hù)及三江源氣候變化與生態(tài)響應(yīng)提供科學(xué)依據(jù)。同時，精確反演地表溫度對于全球氣候變化的研究、地球環(huán)境災(zāi)害監(jiān)測及自然資源管理有著極其重要的意義［2］。

單通道算法是利用衛(wèi)星傳感器上的一個熱紅外通道來獲得地表的輻射能量，并利用大氣的溫度廓線、濕度廓線和壓力廓線，建立大氣輻射傳輸方程，以此為基礎(chǔ)來消除大氣以及地表比輻射率對地表溫度的影響，進(jìn)一步通過普朗克函數(shù)反演得到地表溫度［3］。覃志豪等［4］利用Landsat TM 數(shù)據(jù)的第6波段，推出了只需3個參數(shù)反演地表溫度的算法，其具有較高反演精度。毛克彪等［5］首次利用覃志豪提出的基于Landsat TM 數(shù)據(jù)的單窗反演算法，對MODIS 衛(wèi)星的第31、32 波段建立了輻射率與溫度的線性經(jīng)驗關(guān)系式。毛克彪等［6］首次推導(dǎo)出了適用于ASTER 傳感器的單窗算法，該算法更適合于小區(qū)域的溫度空間差異分析。Landsat 衛(wèi)星具有空間分辨率高、觀測時間序列長的優(yōu)勢，已經(jīng)有接近40 a的觀測數(shù)據(jù)。其次，Landsat 是一種多光譜數(shù)據(jù)，光譜信息豐富且存儲數(shù)據(jù)多，其熱紅外波段對地物的熱信息比較敏感。因此，本文選用Landsat數(shù)據(jù)反演黃河源區(qū)夏季的地表溫度。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

黃河源區(qū)（33°～35°N，95°～102°E）位于青藏高原東部、青海省南部（圖1），南北長約171 km，東西長約232 km，面積約14.75 km2，平均海拔在4000 m左右，最高約6200 m。研究區(qū)屬高原大陸性氣候，氣溫相對較低，早晚溫差大，夏季多年平均氣溫為3.67 ℃。源區(qū)蒸散發(fā)相對強(qiáng)烈，年均蒸發(fā)量達(dá)到1200～2000 mm，2001—2009年黃河源年蒸散量增加明顯［7-8］。

圖1 研究區(qū)示意圖Fig.1 Schematic diagram of the study area

1.2 研究數(shù)據(jù)

本研究選取1990、1995、2000、2005、2010 年7月或8 月Landsat5 的TM 數(shù)據(jù)，選取2015、2018 年7月或8 月Landsat8 的OLI/TIRS 數(shù)據(jù)，由于Landsat 掃描面積相對較小，同一月份很難拼出完整覆蓋研究區(qū)的影像，需要多幅影像鑲嵌而成，選取時云量低于10%，則向相鄰（前后）年度取缺少部分的影像。本文采用空間分辨率為30 m 的SRTM-DEM 高程數(shù)據(jù)，來自中科院地理空間數(shù)據(jù)云網(wǎng)站（http://www.gscloud.cn）；氣象數(shù)據(jù)來自中國氣象數(shù)據(jù)網(wǎng)（http://data.cma.cn）；土壤濕度數(shù)據(jù)來自GLDAS 全球陸面數(shù)據(jù)同化系統(tǒng)（https://disc.gsfc.nasa.gov/datasets）。

1.3 研究方法

1.3.1 反演方法覃志豪針對Landsat 5數(shù)據(jù)的熱紅外波段特點(diǎn)，根據(jù)地表熱輻射傳導(dǎo)方程推導(dǎo)出了精度較高的地表溫度反演算法［4］，主要是將大氣和地表的作用直接反映在公式中，公式如下［4］：

式中：Ts為實(shí)際地表溫度（K）；a、b 為回歸系數(shù)（a=-37.3556，b=0.4586）；Tsensor為傳感器亮度溫度（K）；Ta為大氣平均作用溫度（K）；τ為大氣透過率；ε為地表比輻射率；C、D為關(guān)于τ與ε的函數(shù)。該方法適用于只有一個熱紅外波段的遙感數(shù)據(jù)［9］，且僅需地表比輻射率、大氣透過率、大氣平均作用溫度3個參數(shù)即可實(shí)現(xiàn)地表溫度的反演。

國內(nèi)外學(xué)者針對Landsat8 TIRS的第10波段，對單窗算法公式中的a、b 2 個參數(shù)進(jìn)行擬合，得到Landsat8 TIRS波段10的地表溫度反演算法公式［4］：

式中：Ts10為Landsat8 TIRS 波段10 反演得到的地表溫度（K）；Tsensor10為Landsat8 TIRS 波段10 的傳感器亮度溫度（K）；a10、b10為回歸系數(shù)（表1）；C10、D10為關(guān)于中間變量τ10與ε10的函數(shù)，根據(jù)以下公式得到：

表1 Landsat8 TIRS波段10的回歸系數(shù)取值［10］Tab.1 Regression coefficient value of Landsat8 TIRS band 10［10］

1.3.2 分析方法本文分析利用一元線性回歸趨勢分析方法來計算每個像元在1990—2018年（間隔3～5 a）的變化趨勢，分析黃河源區(qū)夏季地表溫度的變化情況，計算公式為［11］：

式中：QinLSTslope為每個像元在這7 a 變化趨勢的斜率；i為年序號（i=1, 2, 3, …,n）；QinLSTi為第i年夏季的覃志豪算法地表溫度值。QinLSTslope＞0 說明地表溫度在上升，QinLSTslope＜0，說明地表溫度值在下降。

數(shù)學(xué)家皮爾遜基于極大似然法提出變量之間的相關(guān)系數(shù)（r），被用來衡量2 個變量之間線性關(guān)系的強(qiáng)弱。r的正負(fù)代表著變量間的相關(guān)性質(zhì)，若其r越接近于1，則說明兩組數(shù)據(jù)之間的正相關(guān)性越強(qiáng)。本文通過計算皮爾遜相關(guān)系數(shù)來反映地表溫度與其影響因素之間的關(guān)系，計算公式如下［12］：

式中：i為第幾年；n為時間序列長度；xi、yi為第i年的樣本值；xˉ、yˉ分別為樣本值x、y的平均值。

1.3.3 反演參數(shù)的確定覃志豪地表溫度單窗反演算法需要確定3 個物理參數(shù)：（1）地表比輻射率，指物體與黑體在同溫度、同波長下輻射出射度的比值［13-14］。在確定傳感器測距范圍和像素尺寸時，地表比發(fā)射率的大小主要取決于地表材料的組成、形狀和結(jié)構(gòu)［15］。研究區(qū)的地表比輻射率根據(jù)其像元中的土地利用類型及植被指數(shù)計算［16-19］，本研究區(qū)將地表劃分為水體、植被、自然地表（裸地+建筑用地+冰川）。（2）大氣透過率，指電磁波在大氣中傳播時，經(jīng)大氣衰減后的電磁輻射通量與入射時電磁輻射通量的比值，其大小主要由大氣狀況、入射點(diǎn)高度等因素決定。本文利用覃志豪等提出結(jié)合水汽含量動態(tài)變化的大氣透過率方程進(jìn)行估算［20-25］。（3）大氣平均作用溫度，主要受大氣分子的活躍狀態(tài)、大氣剖面溫度及大氣成分的運(yùn)動特征等影響［26］。本文將根據(jù)大氣平均氣溫與地表高程之間的關(guān)系，通過大氣模擬，得出在不同大氣模式及不同大氣剖面的背景下，估算大氣平均作用溫度［27-28］。其中地表比輻射率的計算需要用Landsat 多光譜數(shù)據(jù)求得植被指數(shù)NDVI，再利用像元二分模型計算得到植被覆蓋率，結(jié)合不同地物類型的地表求得整個源區(qū)地表比輻射率。大氣透過率主要受大氣水汽含量的影響，借助經(jīng)驗公式求得大氣透過率的值。根據(jù)不同大氣模式的經(jīng)驗公式估算得到大氣平均作用溫度。

1.3.4 地表溫度反演原始數(shù)據(jù)經(jīng)過一系列預(yù)處理后，需要進(jìn)行輻射亮溫的計算，利用Planck 公式進(jìn)行亮度溫度值的轉(zhuǎn)化，即可獲得地表溫度，公式如下［4］：

式中：Tsensor為傳感器亮度溫度（K）；Lsensor為經(jīng)過預(yù)處理之后的輻射亮度（w·m-2·sr-1·μm-1）；K1、K2為常量，可從影像數(shù)據(jù)文件中獲得（表2）。

表2 輻射常數(shù)KTab.2 Radiation constant K

根據(jù)公式（1）、（4）、（6）即可反演得到地表溫度。

2 結(jié)果與分析

2.1 反演結(jié)果的驗證

將覃志豪反演地表溫度（簡稱反演地表溫度）與MODIS 地表溫度進(jìn)行對比驗證。由于MODIS 數(shù)據(jù)僅在2000 年以后才有數(shù)據(jù)產(chǎn)品，因此選取2000、2005、2010、2015、2018年夏季MODIS 地表溫度數(shù)據(jù)進(jìn)行驗證（圖2）。

圖2 反演地表溫度與MODIS地表溫度對比Fig.2 Comparison between retrieved surface temperature and MODIS surface temperature

總體上，反演地表溫度的平均值高于MODIS 地表溫度的平均值，這是由于相對于其他熱紅外傳感器，Landsat 數(shù)據(jù)具有非常高的空間分辨率，而MODIS 數(shù)據(jù)分辨率較低，其混合像元包含的地物信息變化相對更大，會拉低整體的地表溫度值，另外黃河源區(qū)扎陵湖及鄂陵湖的面積較大，進(jìn)一步影響了黃河源區(qū)夏季地表溫度的平均值。5 個年份的地表溫度平均值在2015 年發(fā)生突變，MODIS 地表溫度平均值在2015 年之前是先下降而后再上升，而反演地表溫度在2015 年之前呈現(xiàn)先上升再下降而后再上升的趨勢，在2015 年之后，兩種數(shù)據(jù)地表溫度的平均值呈下降趨勢。在年際尺度上，兩種地表溫度平均值有所差異，因此有必要對兩種地表溫度數(shù)據(jù)進(jìn)行相關(guān)性分析。

2000—2018 年反演地表溫度與MODIS 地表溫度呈正相關(guān)，其中2015年兩種地表溫度數(shù)據(jù)的相關(guān)性最高，r為0.59，2018 年兩種地表溫度的相關(guān)性最低，其余各年份r值不是很高，在0.43～0.53 之間（圖3），主要原因為兩種地表溫度數(shù)據(jù)的空間分辨率不一致；反演的Landsat數(shù)據(jù)是由不同時間的影響鑲嵌而成，存在一定誤差；少量云的存在會影響反演的結(jié)果。

圖3 反演地表溫度與MODIS地表溫度的相關(guān)性Fig.3 Correlation between retrieved surface temperature and MODIS surface temperature

2.2 地表溫度的時空分布特征

利用覃志豪單窗反演算法反演得到黃河源區(qū)夏季地表溫度空間分布（圖4a～g）。整體上，黃河源區(qū)北部地表溫度最高，可達(dá)50 ℃，這是因為北部地區(qū)海拔較低，且耕地及建筑用地較多；黃河源中部地區(qū)地表溫度呈現(xiàn)藍(lán)色的區(qū)域為冰川，其溫度最低，約為-24 ℃，其他藍(lán)色的區(qū)域有部分是云層覆蓋，導(dǎo)致地表溫度較低呈現(xiàn)藍(lán)色。

圖4 反演地表溫度的空間分布及年際變化趨勢Fig.4 Spatial distribution and interannual variation trend of retrieved surface temperature

局部來看，黃河源區(qū)西部的兩大湖泊，扎陵湖及鄂陵湖的表面溫度相對較低；1995 年黃河源區(qū)南部、2005年及2010年黃河源區(qū)西南部地表溫度有負(fù)值出現(xiàn)（圖4b、d～e），由于零散的云覆蓋所導(dǎo)致。

根據(jù)線性回歸公式（6）計算出7 個年份地表溫度的變化斜率值（圖4h），對地表溫度年平均值序列變化進(jìn)行趨勢分析［29］，從圖中可以看出，1990—2018 年，黃河源區(qū)北部及中部邊緣山區(qū)地表溫度呈現(xiàn)上升趨勢，黃河源西部及東部地區(qū)地表溫度呈現(xiàn)下降趨勢。

2.3 地表溫度影響因素分析

地表溫度主要受高程、土壤濕度、氣候、植被狀況、用地類型以及人類活動等因素的共同影響。黃河源區(qū)建筑用地少，人類活動范圍較小，因此黃河源區(qū)地表溫度主要受自然因素影響。降水和氣溫通過水熱交換影響了每個像元地表溫度的大??；在眾多影響因素中，地表高程與地表溫度的相關(guān)性最高；土壤濕度也影響了地表溫度的空間分布，植被及土地利用類型對地表溫度的影響主要是通過地表比輻射率體現(xiàn)的［30］。因此，探討高程、土壤濕度、氣候、植被狀況和土地類型等因素對黃河源區(qū)夏季地表溫度時空分布的影響十分必要。

2.3.1 高程根據(jù)研究區(qū)高程數(shù)據(jù)，選取高程起伏較大、緯度變化小、地物類型變化不大的區(qū)域，以保證在高程與地表溫度的相關(guān)性分析中，高程是最主要的影響因素。提取其7個年份的高程以及地表溫度并計算其r值（圖5），各年份地表溫度與高程存在明顯的負(fù)相關(guān)關(guān)系，r值平均為-0.65，2000、2015、2018 年均達(dá)0.70 以上（圖5c、f～g）。高程平均每增高100 m，地表溫度下降約0.925 ℃。

圖5 地表溫度與高程的相關(guān)性Fig.5 Correlation between surface temperature and elevation

2.3.2 土壤濕度選用GLDAS 土壤濕度數(shù)據(jù)產(chǎn)品（https://disc.gsfc.nasa.gov），該數(shù)據(jù)的空間分辨率較為粗糙，僅為0.25°，無法在空間上反映地表溫度與土壤濕度的相關(guān)關(guān)系。但是黃河源區(qū)面積大，可以在空間上進(jìn)行對比，可在一定程度上反映土壤濕度與地表溫度的變化情況。黃河源區(qū)土壤濕度的高值區(qū)分布在研究區(qū)的南部及東南部，而低值區(qū)域分布在東北部及西部（圖6a）?？臻g上，黃河源區(qū)東部及西部地表溫度較高，其土壤濕度低；黃河源區(qū)東南部地表溫度相對較低，其土壤濕度高，某種程度上可以說明地表溫度與土壤濕度在空間分布上存在負(fù)相關(guān)關(guān)系（圖6）。

圖6 2015年土壤濕度與地表溫度的空間分布Fig.6 Spatial distributions of soil moisture and surface temperature in 2015

2.3.3 氣溫將氣溫柵格數(shù)據(jù)重新降尺度至30 m空間分辨率，并投影至與地表溫度空間分布圖一致的坐標(biāo)系，計算出黃河源區(qū)夏季地表溫度與夏季氣溫的相關(guān)性空間分布（圖7），整個黃河源區(qū)地表溫度與氣溫呈正相關(guān)的部分約占55%，r平均值約為0.38，主要分布在黃河源區(qū)的西部、北部和東南部山區(qū)。源區(qū)西部及東部存在明顯的分帶現(xiàn)象，是由于不同時相的影像鑲嵌而成。

圖7 地表溫度與氣溫的相關(guān)性Fig.7 Correlation between surface temperature and air temperature

2.3.4 降水量利用ENVI軟件計算得出黃河源區(qū)夏季地表溫度與夏季降水量的相關(guān)性空間分布（圖8），進(jìn)一步統(tǒng)計柵格，其中約70%的像元其r為負(fù)值，其中r小于-0.4 的像元約占20%，說明隨著該地區(qū)降水量的增加，相應(yīng)的黃河源區(qū)地表溫度降低，呈負(fù)相關(guān)關(guān)系的像元主要分布在黃河源區(qū)西部兩大湖泊附近，及黃河源中部山區(qū)、南部東部邊緣地區(qū)。

圖8 地表溫度與降水量的相關(guān)性Fig.8 Correlation between surface temperature and precipitation

2.3.5 用地類型提取不同地物類型的地表溫度（圖9），不難看出，1990—2018年水體的地表溫度基本不變，其他地物的地表溫度均呈現(xiàn)上升的趨勢；冰川與裸地的地表溫度升高最快，冰川從1990年平均地表溫度-7.0 ℃上升到2018 年平均地表溫度3.0 ℃，增幅達(dá)10.0 ℃；裸地由1990 年平均地表溫度23.7 ℃上升到2018 年的34.5 ℃，增幅達(dá)10.8 ℃；耕地、植被的地表溫度呈現(xiàn)上升趨勢，但增速較緩，分別由1990 年平均地表溫度27.2 ℃、21.2 ℃上升到2018年的34.1 ℃、27.3 ℃，增幅分別達(dá)6.9 ℃、6.1 ℃。

圖9 用地類型與地表溫度年際變化Fig.9 Interannual variation of land use type and surface temperature

3 結(jié)論

本文選用7 a 的Landsat 衛(wèi)星數(shù)據(jù)，利用覃志豪單窗算法反演黃河源區(qū)1990—2018 年（間隔3～5 a）夏季的地表溫度，通過一元線性回歸及相關(guān)性分析方法，分析得到黃河源區(qū)地表溫度的時空分布特征和變化趨勢及影響因素，結(jié)論如下：

（1）1990—2018年研究區(qū)北部及中部邊緣山區(qū)地表溫度呈現(xiàn)上升趨勢，研究區(qū)西部及東部地區(qū)地表溫度呈現(xiàn)下降趨勢。

（2）1990—2018 年水體的地表溫度基本不變，其他地物的地表溫度均呈現(xiàn)上升的趨勢；冰川與裸土的地表溫度升高最快，增幅分別達(dá)10.0 ℃、10.8 ℃；耕地及植被的地表溫度值增長速度較為緩慢，由于黃河源區(qū)植被以草地為主，地表溫度對植被覆蓋變化的響應(yīng)存在一定的正相關(guān)關(guān)系。

（3）研究區(qū)高程與地表溫度r為-0.65，即高程平均每升高100 m，地表溫度約下降0.925 ℃；黃河源區(qū)東部及西部地表溫度較高，相應(yīng)的其土壤濕度較低；黃河源區(qū)東南部地表溫度相對較低，其土壤濕度較高。

（4）氣溫的分布對地表溫度的時空分布也有一定影響；降水量與地表溫度呈負(fù)相關(guān)關(guān)系，黃河源區(qū)約70%的像元其r為負(fù)值，說明隨著該地區(qū)降水量的增加，相應(yīng)區(qū)域的地表溫度降低，呈負(fù)相關(guān)關(guān)系的像元主要分布在黃河源區(qū)西部兩大湖泊附近，及黃河源中部山區(qū)、南部東部邊緣地區(qū)。