黃河源永曲河流域地表溫度與土地覆蓋的時空變化特征分析

丁圓圓， 趙健赟*, 楊靜, 趙沁浩, 王祖順, 趙利江

(1.青海大學(xué)地質(zhì)工程系， 西寧 810016； 2.青海省剛察縣氣象站， 剛察 812300； 3.青海省基礎(chǔ)測繪院, 西寧 810001)

地表溫度是地表與大氣相互作用過程中最重要的物理學(xué)參量之一[1]，在全球和局部氣候變化、生態(tài)環(huán)境、城市熱島效應(yīng)等研究中具有重要意義[2-3]，其時空分布與生態(tài)環(huán)境、地形因子、土地覆蓋類型等緊密相關(guān)，因此，關(guān)于地表溫度的反演及其影響因素的研究受到了越來越多學(xué)者的關(guān)注,已經(jīng)成為全球變化領(lǐng)域的研究熱點之一。

目前中外已經(jīng)開展了較多地表溫度領(lǐng)域的研究，大范圍地表溫度數(shù)據(jù)一般通過遙感影像進(jìn)行反演，方法主要有輻射傳輸方程法[4-7]、單窗算法[8-9]、劈窗算法[10-11]。對地表溫度變化及其影響因素方面的研究也已取得了一定的成果，例如，王佳等[12]分別利用普通最小二乘法(ordinary least square, OLS)線性回歸、地理加權(quán)回歸(geographically weighted regression, GWR)模型擬合土地覆蓋比例和地表溫度的關(guān)系，結(jié)果表明不同土地覆蓋類型的地表溫度存在顯著差異，且兩者之間的定量關(guān)系受到空間異質(zhì)性的影響。毛洋等[13]和李軍等[14]認(rèn)為地表溫度的時空分布與地表類型密切相關(guān)。另外，有不少學(xué)者對地表溫度隨地形因子的時空變化規(guī)律也進(jìn)行了研究，趙偉等[15]得到了地表溫度隨海拔高度、坡度的增加而逐漸增加的結(jié)論。

李春波等[16]運用統(tǒng)計學(xué)方法分析阜平縣地表溫度與坡向之間的關(guān)系，結(jié)果表明不同坡向的地表溫度分布具有一定的差異,陽坡的地表溫度高于陰坡。關(guān)于歸一化植被指數(shù)(normalized difference vegetation index, NDVI)與地表溫度的關(guān)系也引起了不少學(xué)者的關(guān)注，但是得到的結(jié)論不完全一致，許劍輝等[17]通過GWR模型擬合二者之間的相關(guān)關(guān)系，得到了NDVI聚集密度對地表溫度有著從正到負(fù)的影響的結(jié)論[17]；Mishra等[18]和Malik等[19]認(rèn)為地表溫度隨植被覆蓋度的變化而變化。張曉娟等[20]則認(rèn)為地表溫度與NDVI呈負(fù)相關(guān)。在青藏高原氣候暖濕化的背景下，黃河源作為高寒生態(tài)系統(tǒng)敏感和脆弱地帶，其地表覆蓋、溫度等環(huán)境氣候因子的變化引起了廣泛的關(guān)注，而以上研究主要集中于低海拔地區(qū)地表溫度的相關(guān)研究中，缺少黃河源高寒草甸覆蓋區(qū)域地表溫度時空變化特征與影響因素的研究。本研究以黃河源區(qū)永曲河流域為例，同時兼顧該區(qū)域植被生長特征、時間尺度、遙感數(shù)據(jù)源的獲取與質(zhì)量等問題，選用2008年、2020年兩期Landat遙感影像反演地表溫度，并與土地覆蓋、地形因子和NDVI的關(guān)系進(jìn)行分析和探討，以獲得研究區(qū)地表溫度的時空變化特征及其規(guī)律。

1 研究區(qū)概況

研究區(qū)為青海省東南部黃河源地區(qū)的永曲河流域，處于101°20′～101°40′E，34°10′～34°30′N，如圖1所示，位于青海省東南部，隸屬黃南藏族自治州，是黃河流域重要的水源地，總面積約為745.76 km2。海拔在3 317～4 451 m，總體地勢表現(xiàn)為東北高西南低，且東部小部分地區(qū)常年有凍土積雪分布。流域氣候為高原大陸性氣候，屬于高原亞寒帶半濕潤氣候區(qū)，由于海拔較高，地勢復(fù)雜和受季風(fēng)影響，高原大陸性氣候特點比較明顯。每年5—10月份溫暖、多雨，11月—次年4月份寒冷、干燥、多大風(fēng)天氣。常年風(fēng)向為西北風(fēng)，最大風(fēng)速可達(dá)到23.7 m/s，年平均風(fēng)速約為2.2 m/s，年平均降水量在400 mm以上，植被類型以草地為主，且僅在高海拔地區(qū)零星分布小灌木叢。

圖1 研究區(qū)遙感影像圖和高程Fig.1 Remote sensing image and altitude in the study area

2 數(shù)據(jù)與研究方法

2.1 數(shù)據(jù)來源與預(yù)處理

本文中使用的數(shù)據(jù)主要包括 Landsat、數(shù)字高程模型(digital elevation model, DEM)、土地利用類型和氣象站數(shù)據(jù)。其中Landsat影像來源于地理空間數(shù)據(jù)云，軌道號為132/36，空間分辨率為30 m，TM5和TM8成像時間分別為2008年8月 17日和2020年8月2日;DEM數(shù)據(jù)來源于地理空間數(shù)據(jù)云，經(jīng)緯度為N34°E101°，空間分辨率為30 m；土地覆蓋數(shù)據(jù)來源于中國科學(xué)院資源環(huán)境科學(xué)與數(shù)據(jù)中心，空間分辨率為30 m。

對收集到的數(shù)據(jù)進(jìn)行了去云、輻射定標(biāo)、大氣校正、投影、重分類等處理，獲得了研究區(qū)2008、2020年的30 m分辨率影像數(shù)據(jù)和DEM數(shù)據(jù)。

2.2 LST反演算法與原理

研究采用輻射傳輸法反演地表溫度，該方法可用于任何熱紅外遙感波段，其基礎(chǔ)理論均基于普朗克定律的熱輻射傳輸方程。為提高地表溫度反演精度，中外學(xué)者們開展了大量的研究，例如，為解決地面物體某些性質(zhì)不能確定的問題，將“黑體”引入輻射傳輸理論；進(jìn)而又考慮到灰體在自然界的普遍存在，將地表比輻射率應(yīng)用于輻射傳輸模型，所需參數(shù)易于獲得，計算過程較為簡單，反演精度高，結(jié)果可靠[21-22]。輻射傳輸法先估計大氣對地表熱輻射的影響，然后把這部分大氣影響從衛(wèi)星傳感器所觀測到的熱輻射總量中減去，進(jìn)而得到地表熱輻射強(qiáng)度，再把這一熱輻射強(qiáng)度轉(zhuǎn)化為相應(yīng)的地表溫度。

2.2.1 熱紅外波段輻射亮度的計算

對Landsat熱紅外波段進(jìn)行輻射定標(biāo)處理，得到熱紅外波段的輻射亮度Lλ[23]為

Lλ=GainDN+Bias

(1)

式(1)中：DN表示傳感器記錄的像元亮度值；Gain表示增益，W/(m2·μm·sr)；Bias表示偏移，W/(m2·μm·sr)。

2.2.2 地表比輻射率的計算

(1)歸一化植被指數(shù)的獲取。歸一化植被指數(shù)NDVI是反映植被生長狀態(tài)、地表覆蓋植被狀況的一種遙感指標(biāo)，可以通過近紅外和紅光波段獲得，即

(2)

式(2)中：NIR為遙感多波段圖像中的近紅外波段；R為紅波段；NDVI位于-1～1，負(fù)值表示地面覆蓋為云、水、雪等，0表示有巖石或裸土等，正值表示有植被覆蓋，且隨覆蓋度增大而增大。

(2)植被覆蓋度的計算。將整個影像的地類大致分為水體、植被和建筑，采用混合像元分解法計算植被覆蓋度FV，即

(3)

式(3)中：NDVIV和NDVIS分別表示植被和裸地的NDVI值。NDVIV=0.7，NDVIS=0.05。當(dāng)NDVI>NDVIV時，F(xiàn)V取值為1；NDVI

(3)地表比輻射率的計算。將水體像元的比輻射率賦值為0.995，自然表面和城鎮(zhèn)像元的比輻射率[24]為

ε自然表面=0.962 5+0.061 4FV-0.046 1FV2

(4)

ε城鎮(zhèn)=0.958 9+0.086FV-0.067 1FV2

(5)

2.2.3 相同溫度下黑體輻射亮度的計算

依據(jù)Landsat的成像時間以及中心經(jīng)緯度得到大氣剖面參數(shù)的向上輻射亮度Lup、大氣向下輻射亮度Ldown、大氣在熱紅外波段的透過率τ，計算B(TS)為

B(TS)=[Lλ-Lup-τ(1-ε)Ldown]/(τε)

(6)

式(6)中：B(TS)表示溫度為TS的黑體輻射亮度。

2.2.4 地表溫度的計算

由普朗克公式的反函數(shù)反演獲得真實地表溫度TS為

TS=k2/ln[k1/B(TS)+1]

(7)

式(7)中：k1、k2為是輻射定標(biāo)參數(shù)，隨傳感器的不同，而有不同的取值。對于Landsat 5數(shù)據(jù),k1=607.76 W/(m2·μm·sr),k2=1 260.56k；k是熱力學(xué)溫度，K。對于Landsat 8數(shù)據(jù)，k1=774.89 W(m2·μm·sr)，k2=1 321.08k。

2.3 地理加權(quán)回歸分析

地理加權(quán)回歸模型通過建立空間范圍內(nèi)每個點處的局部回歸方程，來探索研究對象在某一尺度下的空間變化及相關(guān)驅(qū)動因素，并可用于對未來結(jié)果的預(yù)測。它是對普通線性回歸模型的擴(kuò)展，其顯著的特征是考慮到了空間異質(zhì)性，將數(shù)據(jù)的地理位置嵌入到回歸參數(shù)中，因此具有更高的準(zhǔn)確性，其公式為

i=1,2,…，n

(8)

式(8)中：(ui,vi)是第i個樣點的坐標(biāo)；εi是第i個樣點的隨機(jī)誤差；β0是第i個樣點的截距常量，βk(ui,vi)是第i個樣點的第k個回歸參數(shù)，是關(guān)于地理位置的函數(shù)，在估算的過程中采用權(quán)函數(shù)的方法得到。本文中從局部出發(fā)，以地表溫度(land surface temperature, LST)為因變量，NDVI為解釋變量，并采用高斯函數(shù)作為空間權(quán)重矩陣函數(shù)和修正后的Akaike信息量準(zhǔn)則(Akaike information criterion, AIC)確定帶寬，分析空間異質(zhì)性分布的變量之間的線性關(guān)系。

3 研究結(jié)果與分析

3.1 地表溫度時空變化特征

利用輻射傳輸法對永曲河流域2期Landsat影像做地表溫度反演，結(jié)果如圖2所示，并統(tǒng)計反演結(jié)果的最大值、最小值和均值，結(jié)果如表1所示。從表1可知，2020年地表溫度的最大值、最小值和均值都要高于2008年，其中最小值增加了5.147 ℃，變化最大，這表明2020年的地表溫度相對于2008年整體呈顯著上升趨勢。

圖2 2008年、2020年地表溫度反演結(jié)果Fig.2 Inversion results of LST in 2008 and 2020

表1 2008年、2020年地表溫度統(tǒng)計結(jié)果

3.2 地表溫度影響因素的分析

地表溫度受多方面因素的影響，為了分析永曲河流域地表溫度變化的原因，利用DEM和土地利用類型數(shù)據(jù)，從高程、坡度、坡向、土地覆蓋類型4個方面定量分析地表溫和NDVI的時空分布和變化特征。

3.2.1 地表溫度的變化特征

根據(jù)研究區(qū)的地形特征，將各地形因子按以下方式進(jìn)行等級劃分。

(1)將海拔高度每隔400 m分為一個等級，研究區(qū)的海拔高度劃分為3 317～3 736 m、3 736～3 959 m、3 959～4 451 m共3個等級，如圖3所示。

圖3 高程分級圖Fig.3 Image of altitude classification

(2)將坡度劃分為平坡 (0°～5°)、緩坡 (5°～15°)、斜坡(15°～25°)、陡坡(25°～35°)和急坡(35°～90°)共5個等級，如圖4所示。

圖4 坡度分級圖Fig.4 Image of slope classification

(3)坡向以正北方向為0°，將0°～360°劃分為平緩坡(0°)、陰坡(0°～45°和315°～360°)、半陰坡(45°～90°和270°～315°)、陽坡(135°～225°)、半陽坡(90°～135°和225°～270°)共5個等級，如圖5所示。

圖5 坡向分級圖Fig.5 Image of aspect classification

對不同高程、坡度、坡向等級下的地表溫度和NDVI進(jìn)行統(tǒng)計，如表2、表3所示,統(tǒng)計表明：

表2 2008年地表溫度和歸一化植被指數(shù)分類統(tǒng)計Table 2 Classification of LST and NDVI in 2008

表3 2020年地表溫度和歸一化植被指數(shù)分類統(tǒng)計Table 3 Classification of LST and NDVI in 2020

(1)海拔與地表溫度、NDVI呈現(xiàn)負(fù)相關(guān)的趨勢，即海拔越高的地方，地表溫度和NDVI越低；當(dāng)海拔低于3 959 m時，隨著海拔升高，地表溫度和NDVI下降幅度較小，海拔大于3 959 m時，地表溫度和NDVI下降較快，說明二者受海拔影響更顯著。

(2)坡度也會對地表溫度和NDVI造成一定的影響：地表溫度和NDVI隨坡度增加而降低，平坡處的地表溫度和NDVI是最高的，急坡處的地表溫度和NDVI最低，平坡處溫和的環(huán)境有利于植物的生長，植物長勢良好，相應(yīng)地，該處的NDVI也較高；

(3)按照陽坡、平緩坡、半陽坡、半陰坡、陰坡的順序，地表溫度和NDVI逐漸降低；2008年陰坡和陽坡之間的地表溫度相差1.169 ℃，2020年陰坡和陽坡之間的地表溫度相差0.626 ℃，其他坡向之間的差異較小。

3.2.2 不同土地覆蓋類型的地表溫度變化特征

地表溫度與土地覆蓋類型緊密相關(guān)，為分析兩者之間的關(guān)系，將研究區(qū)的土地覆蓋類型與地表溫度做分區(qū)統(tǒng)計,如表4所示。由于研究區(qū)特殊的地理位置和環(huán)境，土地覆蓋主要為草地和小灌木叢兩大類，研究區(qū)內(nèi)低海拔地區(qū)以草地覆蓋為主，周邊高海拔地區(qū)僅有小灌木叢零星分布。統(tǒng)計結(jié)果表明：低溫區(qū)的分布與小灌木叢的分布在空間上基本保持一致；小灌木叢主要分布在3 959～4 451 m的高海拔地區(qū)，草地主要分布在3 736～3 959 m的低海拔地區(qū)，高寒環(huán)境下，小灌木叢的生長和覆蓋度受到抑制，從而造成了小灌木叢的NDVI要比草地的NDVI低的現(xiàn)象，這說明NDVI受海拔影響很大，且二者呈現(xiàn)出負(fù)相關(guān)的關(guān)系。在2008、2020年，小灌木叢處的地表溫度都較低，溫度最大相差2.043 ℃，而草地處的地表溫度變化則很小，但總體上溫度呈現(xiàn)上升的趨勢。

表4 不同土地覆蓋類型的地表溫度和歸一化植被指數(shù)分類統(tǒng)計Table 4 Classification of LST and NDVI for different land cover types

將反演的地表溫度分為低溫區(qū)(<21 ℃)、常溫區(qū)(21～27 ℃)、高溫區(qū)(>27 ℃)，并分析2008年、2020年的不同土地覆蓋類型占各溫度區(qū)間的情況，如表5所示，對比2個年份不同地類的溫度區(qū)間面積，按照低溫區(qū)、常溫區(qū)、高溫區(qū)的順序，2008年各溫度區(qū)間所占研究區(qū)的面積比例分別為9.147%、77.163%、13.69%，2020年各溫度區(qū)間所占研究區(qū)的面積比例分別為3.583%、90.374%、6.043%，從這兩年不同溫度區(qū)間所占的比例來看，2020年常溫區(qū)面積比例增加了13.211%，低溫區(qū)間和高溫區(qū)間下降的幅度基本相當(dāng)，進(jìn)一步表明研究區(qū)地表溫度呈現(xiàn)上升的趨勢，這與2008年、2020年的研究區(qū)地表溫度反演結(jié)果保持一致。

表5 不同土地覆蓋類型在各溫度區(qū)間所占的比例Table 5 Proportion of different land cover types in each temperature range

3.3 地表溫度與NDVI的地理加權(quán)回歸分析

3.3.1 NDVI空間格局及變化特征

NDVI是反映地表植被覆蓋狀況的一種遙感指標(biāo)，包含了研究區(qū)的植被信息、水體信息、建筑信息等。本文中分析了2008—2020年的研究區(qū)NDVI的分布狀況，如圖6所示。由圖6中可知， 2個年份的NDVI均具有明顯的空間分布特征，即大部分地區(qū)的NDVI值都比較高，只有東部地區(qū)的NDVI較低；從2個年份的NDVI統(tǒng)計值來看，NDVI呈現(xiàn)出上升的趨勢，即2020年研究區(qū)植被覆蓋度相對于2008年有所增加，且在研究區(qū)的東北部和西南部增加顯著。

圖6 2008年、2020年歸一化植被指數(shù)空間分布Fig.6 The space distribution of NDVI in 2008 and 2020

3.3.2 NDVI對氣溫的影響特征

對2008年和2020年的LST與NDVI構(gòu)建GWR模型，并對回歸系數(shù)做地圖可視化,如圖7所示。從圖7中可以看出：NDVI回歸系數(shù)為正值，這表示NDVI對地表溫度為正影響；2008年NDVI回歸系數(shù)在0.37～9.77變化，2020年的NDVI回歸系數(shù)主要在7.28～8.62，這說明2008年NDVI對地表溫度的影響程度的變化較大，2020年NDVI對LST的影響程度比較集中，變化較??；從回歸系數(shù)的空間變化來看，整體上回歸系數(shù)從西到東逐漸增加，即NDVI對地表溫度的影響逐漸增強(qiáng)。

圖7 2008年、2020年歸一化植被指數(shù)的回歸系數(shù)Fig.7 Regression coefficients of NDVI in 2008 and 2020

擬合優(yōu)度(LocalR2)表示自變量對因變量的解釋程度，研究獲得的LocalR2的空間分布如圖8所示,較大值主要分布在東部地區(qū)，較小值分布在西部地區(qū)，且2個年份的LocalR2值均小于0.5，這表示NDVI對LST的解釋程度較低，受空間異質(zhì)性影響較強(qiáng)烈，也即是除了NDVI影響地表溫度外，還有其他因素對地表溫度有影響。

圖8 2008年、2020年地理加權(quán)回歸模型局部擬合度Fig.8 Degree of local fit of GWR models in 2008 and 2020

3.3.3 GWR殘差分析

以1 km×1 km尺度為例，以LST為因變量，NDVI為解釋變量，對2008年、2020年的回歸方程殘差進(jìn)行分析，結(jié)果如圖9所示，可以看出，研究區(qū)東南部地區(qū)的標(biāo)準(zhǔn)化殘差超過2.5倍，擬合效果不太理想?？紤]到殘差的空間自相關(guān)性會對分析結(jié)果產(chǎn)生一定的影響，為進(jìn)一步分析和驗證結(jié)果，使用全局Moran’sI(莫蘭指數(shù))檢驗殘差的空間自相關(guān)性，經(jīng)計算，2008年和2020年的Moran’sI分別為0.536和0.492，這說明LST與NDVI的關(guān)系受空間異質(zhì)性影響十分強(qiáng)烈。

圖9 2008年、2020年地理加權(quán)回歸模型殘差分布Fig.9 Distribution of GWR residuals in 2008 and 2020

另外，目前關(guān)于地表溫度和NDVI的研究主要通過相關(guān)系數(shù)和普通回歸模型(如最小二乘法)，但在研究區(qū)不同的空間位置， NDVI的變化對地表溫度的影響程度可能并不一致。本文研究使用Moran’sI和地理加權(quán)回歸分別從全局和局部的角度來探討地表溫度與NDVI的相關(guān)關(guān)系。由于地理加權(quán)回歸引入了空間權(quán)重，能有效地定量研究空間非平穩(wěn)性和較好地解決研究對象的空間異質(zhì)性問題，避免了傳統(tǒng)線性回歸模型不能反映回歸參數(shù)的真實空間特征的局限性。

4 結(jié)論

利用輻射傳輸法，反演了永曲河流域2008年、2020年的地表溫度，并分析了地表溫度的時空變化特征及其與地形因子、土地覆蓋類型和NDVI的相關(guān)關(guān)系，得到了以下結(jié)論。

(1)2020年研究區(qū)地表溫度相對于2008年總體上存在較為顯著的上升趨勢。

(2)地表溫度和NDVI受地形因子影響表現(xiàn)出不同的空間分布特征：平均溫度和NDVI隨著海拔和坡度的增加而逐漸降低，且海拔越高負(fù)相關(guān)性越顯著；不同坡向處的地表溫度和NDVI也存在一定的差異，且二者在陽坡和陰坡處的差異更顯著。

(3)研究區(qū)2008年、2020年NDVI與地表溫度呈正相關(guān)關(guān)系。地理加權(quán)回歸殘差具有較強(qiáng)的空間自相關(guān)性，NDVI與地表溫度的關(guān)系受空間異質(zhì)性影響較強(qiáng)烈，而地理加權(quán)回歸法充分考慮了對象的權(quán)重問題，能有效獲得對象的非平穩(wěn)關(guān)系和異質(zhì)性特征，具有一定的優(yōu)勢。

研究結(jié)果對黃河源的生物多樣性、植被保護(hù)、水資源安全問題等有重要意義，也能為改善該區(qū)脆弱的生態(tài)環(huán)境提供一定的參考。本文中采用輻射傳輸法反演地表溫度，仍缺乏地面驗證，且GWR模型中只分析了1 km×1 km尺度上NDVI與地表溫度的關(guān)系，對于二者在不同空間尺度上的定性定量分析有待進(jìn)一步研究。