李建明，馬燕飛，李仁杰，2，4，商國營，李明明

(1.河北師范大學(xué) 資源與環(huán)境科學(xué)學(xué)院，石家莊 050024；2.河北省環(huán)境變化遙感識別技術(shù)創(chuàng)新中心，石家莊 050024；3.邯鄲學(xué)院 地理系，河北 邯鄲 056005；4.河北省環(huán)境演變與生態(tài)建設(shè)實(shí)驗(yàn)室，石家莊 050024)

0 引言

陸地表面溫度(land surface temperature，LST)是陸地表面水熱平衡的重要分量，直接驅(qū)動地表與大氣之間能量交換。觀測與反演LST，對全球氣候變化評估、生態(tài)問題控制、水熱循環(huán)的進(jìn)一步研究具有重要意義。同時(shí)，對農(nóng)業(yè)指導(dǎo)、旱情監(jiān)測、提高農(nóng)田水資源利用率做出了重要貢獻(xiàn)[1]。目前，LST已經(jīng)廣泛應(yīng)用到多種領(lǐng)域，如城市熱島[2]、土壤水估算[3]、蒸散發(fā)(evapotranspiration，ET)估計(jì)[4]等。但現(xiàn)有的地表溫度產(chǎn)品多為大網(wǎng)格尺度，難以滿足田塊尺度的地表熱通量研究[5]。MODIS熱紅外遙感地表溫度產(chǎn)品可以大面積獲取地表熱通量數(shù)據(jù)，但是空間分辨率較低，在中緯度地區(qū)受云雨天氣影響，熱紅外地表溫度空間上不連續(xù)、缺失大量像元，實(shí)際應(yīng)用價(jià)值低[6-7]。近年來，星載熱紅外傳感器空間分辨率與時(shí)間分辨率不能共存的缺陷日趨明顯。在地表溫度反演領(lǐng)域引入LST降尺度方法，使獲取不同時(shí)空分辨率的地表溫度產(chǎn)品成為可能[8]，因此，降尺度技術(shù)就成為連接大尺度LST與區(qū)域小尺度地表熱通量的紐帶。

地表溫度空間降尺度是在低空間分辨率像元值的基礎(chǔ)上，融合高空間分辨率的解釋向量，獲取更多的空間細(xì)節(jié)信息。降尺度根據(jù)算法原理的不同可以分為兩類。一類是基于物理機(jī)制的圖像融合降尺度方法，包括空間濾波分析法[9]、相關(guān)系數(shù)法(相關(guān)統(tǒng)計(jì)分析)[10]、主成分分析法等[11]。二類是基于尺度因子統(tǒng)計(jì)關(guān)系的降尺度方法[12]。Kustas等[13]建立陸地表面溫度與植被指數(shù)的統(tǒng)計(jì)回歸關(guān)系，融合高分辨率的航空數(shù)據(jù)(24 m)進(jìn)行地球同步運(yùn)行環(huán)境衛(wèi)星數(shù)據(jù)降尺度，得到田塊尺度的空間細(xì)節(jié)信息，并提出使用熱紅外遙感與NDVI進(jìn)行TsHARP算法融合，對地表溫度反演；Ma等[14]使用NDVI在Kustas的TsHARP算法基礎(chǔ)上進(jìn)行改進(jìn)，對地表溫度進(jìn)行降尺度。但是，傳統(tǒng)算法難以表征地表溫度與地表參數(shù)之間復(fù)雜統(tǒng)計(jì)關(guān)系。機(jī)器學(xué)習(xí)在地表溫度數(shù)據(jù)降尺度的探索、分析、預(yù)測和模型的建立與評估方面優(yōu)勢明顯，可以在眾多解釋向量中提取需要的信息，獲取高分辨率地表溫度影像，并重建熱紅外地表溫度的缺失像元。孟楊繁宇[15]使用機(jī)器學(xué)習(xí)的方法對黑河流域的蒸散發(fā)進(jìn)行回歸、預(yù)測。受到啟發(fā)，毛克彪等[16]嘗試將機(jī)器學(xué)習(xí)引入LST降尺度研究，使用神經(jīng)網(wǎng)絡(luò)算法對地表溫度進(jìn)行反演，引入AMSR-E被動微波溫度數(shù)據(jù)，并重建熱紅外數(shù)據(jù)在云雨天氣的數(shù)據(jù)缺失。Hutengs等[17]使用隨機(jī)森林算法對1 000 m分辨率的地表溫度進(jìn)行降尺度，利用紅光、近紅外波段的地表反射率和數(shù)字高程模型DEM構(gòu)建最優(yōu)解釋向量集，建立簡單隨機(jī)森林降尺度模型(basic-RF)，得到250 m分辨率的地表溫度。

本文以海河流域?yàn)檠芯繀^(qū)，MODIS地表溫度產(chǎn)品為數(shù)據(jù)源，通過隨機(jī)森林算法對地表溫度進(jìn)行降尺度，建立地表溫度與解釋變量(植被指數(shù)、空氣溫度、下行短波輻射)的降尺度模型。建立了精細(xì)圖像與粗圖像之間的非線性映射關(guān)系，選取植物生長月份進(jìn)行地表溫度的回歸預(yù)測。探討隨機(jī)森林降尺度算法在不同下墊面的表現(xiàn)，刻畫田塊尺度的地表溫度在海河流域的時(shí)空分布格局，為海河流域地表水熱循環(huán)研究提供數(shù)據(jù)參考。

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

1.1 研究區(qū)概況

海河流域位于112°E～120°E、35°N～43°N之間。東臨渤海灣，西倚太行山，南臨黃河，北接蒙古高原[18]。行政區(qū)劃包括京津冀大部、山西省東部、河南省北部及山東省東北部地區(qū)。流域總面積32.06×104km2，占全國總面積的3.3%。海河流域包括五大支流(潮白河、永定河、大清河、子牙河、南運(yùn)河)和一個(gè)小支流(北運(yùn)河)。海河流域地勢西高東低，以高原、山地以及平原等三種地貌為主。海河流域地處溫帶半濕潤半干旱大陸性季風(fēng)氣候區(qū)，流域年平均氣溫在1.5～14.0 ℃，年平均相對濕度在50%～70%之間，年平均降水量約為535 mm，年平均陸面蒸發(fā)量約為470 mm，水面蒸發(fā)量約為1 100 mm。近年來，海河流域氣溫呈明顯上升趨勢，水資源呈減少趨勢[19]。

1.2 數(shù)據(jù)來源與處理

本文主要使用MODIS產(chǎn)品中2～3級的V005版標(biāo)準(zhǔn)數(shù)據(jù)產(chǎn)品，包括地表溫度產(chǎn)品、植被指數(shù)NDVI產(chǎn)品，具體信息見表1。其中，地表溫度產(chǎn)品MOD11A1數(shù)據(jù)集的DN值(有效數(shù)據(jù))范圍為7 500～65 535，地表真實(shí)溫度值由DN值轉(zhuǎn)換得到，轉(zhuǎn)換如式(1)所示。

Ts=(DN)*(scale factor)+offset(K)

(1)

式中：scale factor為比例系數(shù)，取值為0.02;offset為偏移量，取值為0。

植被指數(shù)產(chǎn)品存儲信息與地表溫度產(chǎn)品類似。MODIS數(shù)據(jù)產(chǎn)品是以HDF-EOS格式存儲，正弦格式(sinusoidal)方法進(jìn)行投影的。如圖1所示，本文所需完整海河流域的MODIS數(shù)據(jù)是由四幅影像鑲嵌得到，需利用MODIS數(shù)據(jù)預(yù)處理工具M(jìn)RT(MODIS reprojection tool)對這些產(chǎn)品進(jìn)行重采樣和鑲嵌。首先，考慮到海河流域的地理區(qū)位與空間形狀特征，將影像均采用Albers等面積投影；然后，使用最鄰近法將影像均重采樣為250 m×250 m、500 m×500 m、1 000 m×1 000 m三種柵格大??；最后，使用海河流域矢量邊界將重采樣后的柵格數(shù)據(jù)裁剪到研究區(qū)范圍大小。

圖1 研究區(qū)概況圖、MODIS數(shù)據(jù)分幅及周邊氣象臺站

表1 MODIS數(shù)據(jù)產(chǎn)品

1.3 地面站點(diǎn)數(shù)據(jù)

“海河流域多尺度地表通量與氣象要素觀測數(shù)據(jù)集”[20]是國家青藏高原科學(xué)數(shù)據(jù)中心(https：//data.tpdc.ac.cn/zh-hans/)發(fā)布的2008—2010年自動氣象站觀測數(shù)據(jù)。用戶可通過“寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心”申請并下載該數(shù)據(jù)集。該數(shù)據(jù)集不僅可以定量揭示海河流域主要下墊面地表水熱通量交換特征，同時(shí)也為遙感估算地表蒸散量的驗(yàn)證提供地面“相對真值”。如表2與圖2所示，數(shù)據(jù)來源主要包括三個(gè)自動氣象站點(diǎn)：海河流域密云站(山區(qū)林地)、館陶站(南部農(nóng)田)和大興站(中部城郊)。密云站下墊面是果園(李子樹和蘋果樹)，館陶站下墊面是耕地(玉米、小麥、棉花)，大興站下墊面是耕地(玉米、小麥)，氣象站點(diǎn)下墊面性質(zhì)較為均一。自動氣象站的采集頻率為10 s，且10 min輸出一次。

表2 觀測點(diǎn)數(shù)據(jù)

圖2 氣象站點(diǎn)

1.4 其他數(shù)據(jù)

本文還收集了用于地表溫度降尺度的海河流域部分地區(qū)的空氣溫度、下行短波輻射數(shù)據(jù)。用到的相關(guān)專題圖有：海河流域邊界矢量圖、海河流域土地利用(GlobeLand30 http：//www.globallandcover.com)等。在具體應(yīng)用過程中還對這些專題圖件做了進(jìn)一步的裁剪、重采樣及范圍匹配等處理。

2 研究方法

2.1 隨機(jī)森林原理

隨機(jī)森林是用途廣泛的一種機(jī)器學(xué)習(xí)算法，隨機(jī)森林模型的本質(zhì)是對決策樹算法的改進(jìn)。隨機(jī)森林模型是一個(gè)樹型分類器，輸入一個(gè)二維的矩陣，構(gòu)建沒有剪枝的回歸決策樹，每個(gè)葉子節(jié)點(diǎn)代表一個(gè)回歸結(jié)果，決定了單棵樹的生長過程，隨機(jī)森林的回歸模型采用單棵樹輸出結(jié)果的簡單平均得到。在估算數(shù)據(jù)的過程中，解釋變量的構(gòu)建對于目標(biāo)數(shù)據(jù)相當(dāng)重要。隨機(jī)森林是一個(gè)決策樹的集合，它的功能更多依賴于決策樹的功能[21]，是在其基礎(chǔ)上的擴(kuò)展，從訓(xùn)練庫中提取實(shí)例進(jìn)行判斷，根據(jù)屬性最終決定拿出哪個(gè)樣本來預(yù)測最優(yōu)結(jié)果。因此可以說，隨機(jī)森林是一個(gè)元估計(jì)器，適用于大量子樣本基礎(chǔ)上的決策樹分析，主要使用平均值來提高預(yù)測精度和控制過度擬合[22]，最終的預(yù)測結(jié)果不會擺脫訓(xùn)練樣本的范圍。本文構(gòu)建的模型本質(zhì)上屬于回歸問題，對于回歸問題隨機(jī)森林模型給出所有決策樹預(yù)測結(jié)果的平均。在地表溫度降尺度過程中，地表溫度與各地表參數(shù)(解釋向量)間呈非線性關(guān)系，而隨機(jī)森林降尺度模型對多元共線性不敏感。

本研究采用R語言中random forest數(shù)據(jù)包構(gòu)建隨機(jī)森林降尺度模型，訓(xùn)練樣本為低空間分辨率(500 m、1 000 m)的遙感影像，地表溫度作為因變量，植被指數(shù)、空氣溫度、下行短波輻射作為解釋變量，以獲取更高空間分辨率的地表溫度。

如圖3所示，解釋變量集與地表溫度訓(xùn)練出的某種映射關(guān)系構(gòu)成了隨機(jī)森林地表溫度降尺度模型，即在低空間分辨率影像(1 000 m)上建立地表溫度與地表參數(shù)(解釋變量)的統(tǒng)計(jì)關(guān)系，并假設(shè)這種統(tǒng)計(jì)關(guān)系不隨尺度的大小而變化。計(jì)算如式(2)所示。

LSTpre(EVI)=H(EVI)

(2)

式中：LSTpre表示地表溫度預(yù)測值；EV表示多個(gè)地表參數(shù)組成的最優(yōu)解釋變量集；I表示高空間分辨率影像(250 m)；H表示地表溫度與最優(yōu)解釋變量集的回歸模型。

圖3 降尺度示意圖

2.2 降尺度模型構(gòu)建

隨機(jī)森林地表溫度降尺度由數(shù)據(jù)模型驅(qū)動，其本質(zhì)是以大量解釋向量樣本為基礎(chǔ)，通過隨機(jī)森林模型方法訓(xùn)練樣本，進(jìn)而挖掘遙感圖像更加精細(xì)的空間細(xì)節(jié)。因此，解釋向量集的選擇十分重要。經(jīng)過實(shí)驗(yàn)，本文將表3中的歸一化植被指數(shù)、空氣溫度與下行短波輻射構(gòu)成解釋向量集。本文地表溫度降尺度的核心是：隨機(jī)森林降尺度模型的建立與解釋變量集的選擇。降尺度模型構(gòu)建的具體步驟為：①M(fèi)ODIS LST數(shù)據(jù)與解釋變量數(shù)據(jù)的鑲嵌、重投影等預(yù)處理，并重分類為高分辨率(250 m)和低空間分辨率(500 m、1 000 m)的數(shù)據(jù)集；②將不同分辨率的解釋向量數(shù)據(jù)集、地表溫度數(shù)據(jù)集轉(zhuǎn)換成隨機(jī)森林算法可以識別的CSV文件；③將MODIS LST、解釋向量集作為降尺度模型的訓(xùn)練數(shù)據(jù)進(jìn)行訓(xùn)練；④把高分辨率(250 m)的解釋向量數(shù)據(jù)輸入到訓(xùn)練好的降尺度模型，獲取高分辨率(250 m)的地表溫度預(yù)測結(jié)果；⑤將降尺度模型預(yù)測結(jié)果CSV文件轉(zhuǎn)換回遙感圖像。技術(shù)線路如圖4所示。

表3 解釋變量集類型與名稱

圖4 技術(shù)流程圖

3 降尺度結(jié)果

3.1 降尺度模型篩選

對隨機(jī)森林降尺度模型輸入?yún)?shù)進(jìn)行調(diào)整，構(gòu)建500 m、1 000 m兩個(gè)不同分辨率層級的訓(xùn)練樣本，使用均方根誤差(RMSE)、決定系數(shù)(R2)、偏差(bias)、相對精度(RA)對模型進(jìn)行篩選，確定最優(yōu)模型算法。分別在500 m、1 000 m分辨率層級上使用隨機(jī)森林模型進(jìn)行訓(xùn)練。500 m與1 000 m分辨率層級模型的驗(yàn)證結(jié)果如表4所示。在不同分辨率層級，相同的解釋向量組合的訓(xùn)練集下，隨機(jī)森林降尺度模型的性能表現(xiàn)優(yōu)秀，可以使用更少的數(shù)據(jù)來描述解釋向量集與陸表溫度之間的非線性關(guān)系。在500 m和1 000 m隨機(jī)森林訓(xùn)練模型中，隨機(jī)森林算法在1 000 m分辨率層級下表現(xiàn)出的空間異質(zhì)性更小，訓(xùn)練模型表現(xiàn)更優(yōu)，下一步采用1 000 m分辨率的訓(xùn)練模型進(jìn)行陸表溫度反演。

表4 隨機(jī)森林在500 m、1 000 m層級下模型反演LST的評分

3.2 降尺度模型RF的站點(diǎn)直接驗(yàn)證

地面觀測實(shí)驗(yàn)的目的是為了獲取地面真實(shí)值，以便對遙感估算結(jié)果進(jìn)行驗(yàn)證，但是由于儀器觀測的不確定性和誤差的客觀存在，所以必須對觀測數(shù)據(jù)進(jìn)行處理與質(zhì)量控制[23]?；贛ODIS MOD11A1遙感數(shù)據(jù)來獲取地表溫度，降尺度后空間分辨率為250 m，這與自動氣象儀的觀測尺度是不匹配的。本文地表溫度觀測數(shù)據(jù)采用四分量輻射儀計(jì)算得到的結(jié)果。四分量輻射有效探測范圍分別是：館陶站15.6 m、大興站28 m、密云站30.76 m。一定程度上克服了站點(diǎn)尺度監(jiān)測范圍小的缺陷，計(jì)算如式(3)所示。

(3)

式中：σ為斯蒂芬-玻爾茲曼常數(shù)(5.67×10-8W/m2K4)；ε為地表比輻射率；Rlu為大氣上行長波輻射；Rld為大氣下行長波輻射。

大興站、密云站、館陶站的降尺度LST與輻射四分量對比，平均誤差分別為3.86 K、2.37 K、1.96 K。地面氣象數(shù)據(jù)為點(diǎn)尺度的數(shù)據(jù)，對衛(wèi)星數(shù)據(jù)進(jìn)行驗(yàn)證時(shí)理想的下墊面為勻質(zhì)下墊面。其中，大興站為非勻質(zhì)下墊面(旱地與果園)，因此在尺度匹配時(shí)出現(xiàn)系統(tǒng)誤差，導(dǎo)致降尺度地表溫度與地面氣象數(shù)據(jù)誤差較大；密云站、館陶站為勻質(zhì)下墊面(旱地)，降尺度地表溫度與地面氣象數(shù)據(jù)誤差較小。

圖5 RF 模型的反演的LST與地面站點(diǎn) LST散點(diǎn)圖

3.3 RF降尺度模型在不同下墊面的表現(xiàn)

為了定量描述隨機(jī)森林降尺度模型在不同下墊面的降尺度效果，如圖6(a)所示，本文在海河流域選取了五個(gè)不同下墊面性質(zhì)的子研究區(qū)。其中，子研究區(qū)A下墊面性質(zhì)為草地；子研究區(qū)B下墊面性質(zhì)為林地；子研究區(qū)C下墊面性質(zhì)為旱地；子研究區(qū)D下墊面性質(zhì)為旱地與建設(shè)用地；子研究區(qū)E下墊面性質(zhì)為林草交錯(cuò)地。將五個(gè)子研究區(qū)降尺度結(jié)果升尺度至1 000 m，分別與MODIS原始地表溫度影像對比在不同下墊面性質(zhì)區(qū)域的降尺度效果并分析其空間分布特征。通過表5與圖6(c)可以看出，RF隨機(jī)森林方法在海河流域地表溫度降尺度上總體表現(xiàn)優(yōu)秀，這是由于RF隨機(jī)森林主要使用平均值來提高預(yù)測精度和控制過度擬合，可以更好地捕捉解釋向量與LST之間的函數(shù)關(guān)系，從而獲得高精度的降尺度結(jié)果。在植被覆蓋度高的草地、林地、林草交錯(cuò)地帶降尺度效果R2、RMSE、bias優(yōu)于植被覆蓋度低的建設(shè)用地、旱地，降尺度的相對精度也更高，草地、林地、林草交錯(cuò)地與原始溫度反演結(jié)果一致性更高。與植被覆蓋度高的子研究區(qū)相比旱地、建設(shè)用地RF隨機(jī)森林降尺度優(yōu)勢并不明顯。從影像的RMSE來看，建設(shè)用地達(dá)到3.894 K，說明解釋變量集(NDVI、空氣溫度、下行短波輻射)與LST的關(guān)系可能呈現(xiàn)出不穩(wěn)定的非線性關(guān)系。

從圖6(b)可以看出，在五個(gè)子研究區(qū)內(nèi)，降尺度影像與MODIS原始影像的地表溫度在空間分布上一致，但總體來說降尺度影像的地表溫度要略低于MODIS原始影像。在降尺度地表溫度產(chǎn)品空間分布格局上，能夠合理地反映下墊面溫度變化，對地表溫度的空間分析具有重要意義。

注：A、B、C、D、E代表五個(gè)子研究區(qū)；下標(biāo)1表示降尺度地表溫度；下標(biāo)2表示MODIS原始地表溫度。圖6 RF降尺度模型效果

表5 RF模型反演LST的評分

3.4 重建地表溫度缺失像元

如圖7所示，通過RF|1 000 m|模型的降尺度的地表溫度數(shù)據(jù)，融合250 m分辨率的無缺失像元解釋向量集得到的250 m海河流域地表溫度數(shù)據(jù)，填補(bǔ)了MODIS MOD11A1地表溫度產(chǎn)品由云導(dǎo)致的地表溫度缺失像元。獲取海河流域地表溫度更精細(xì)的空間信息，得到精度更高的地表溫度產(chǎn)品，滿足海河流域地表溫度時(shí)空格局的研究。

注：①代表2008年第241天；②代表2010年第240天。圖7 海河流域降尺度結(jié)果

4 結(jié)束語

本文以海河流域?yàn)檠芯繀^(qū)，探索了利用植被指數(shù)、空氣溫度、下行短波輻射等解釋向量對地表溫度空間降尺度的方法。通過對隨機(jī)森林算法在兩個(gè)空間尺度(500 m、1 000 m)下的訓(xùn)練模型進(jìn)行分析，隨機(jī)森林1 000 m算法表現(xiàn)優(yōu)于隨機(jī)森林500 m。將MODIS MOD11A1地表溫度產(chǎn)品與隨機(jī)森林模型降尺度的地表溫度在五個(gè)子研究區(qū)進(jìn)行驗(yàn)證分析，結(jié)果表明降尺度地表溫度的結(jié)果空間分布較為合理。將降尺度地表溫度與地面氣象站點(diǎn)的輻射四分量溫度進(jìn)行直接驗(yàn)證，在大興站、密云站、館陶站的平均誤差分別為：3.86 K、2.37 K、1.96 K。最終，構(gòu)建了隨機(jī)森林算法地表溫度降尺度模型。將MODIS MOD11A1地表溫度降尺度到250 m分辨率，并利用降尺度地表溫度填補(bǔ)MODIS地表溫度的缺失像元，避免了MODIS MOD11A1地表溫度產(chǎn)品空間分辨率較低、易受天氣影響、有效像元信息缺失嚴(yán)重、實(shí)際應(yīng)用價(jià)值低等問題。海河流域的降尺度地表溫度數(shù)據(jù)為數(shù)據(jù)稀疏區(qū)域的數(shù)據(jù)獲取提供科學(xué)參考。

本文總體上亮點(diǎn)有三點(diǎn)。一是將植被指數(shù)、空氣溫度、下行短波輻射組成的最優(yōu)解釋向量集引入到MODIS地表溫度降尺度，將地表溫度的空間分辨率由1 000 m降尺度到250 m，獲取到更加精細(xì)的地表溫度空間細(xì)節(jié)。二是在海河流域選出五個(gè)子研究區(qū)，分析了在不同下墊面性質(zhì)下，隨機(jī)森林降尺度模型的效果。三是采用輻射四分量溫度對地表溫度降尺度結(jié)果進(jìn)行直接驗(yàn)證，一定程度上改善了驗(yàn)證數(shù)據(jù)與降尺度地表溫度像元尺度不匹配的問題。

但文中仍存在一些不足需進(jìn)一步研究，如通過遙感反演得到的地表溫度具有一定的方向(觀測角度)性，模型沒有考慮衛(wèi)星傳感器傾斜觀測角度造成的地表溫度差異。如何將傾斜觀測的地表溫度歸一到垂向觀測值，是未來研究需進(jìn)一步考慮的內(nèi)容。