基于隨機森林的海河流域地表溫度降尺度

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

(1.河北師范大學 資源與環境科學學院，石家莊 050024；2.河北省環境變化遙感識別技術創新中心，石家莊 050024；3.邯鄲學院 地理系，河北 邯鄲 056005；4.河北省環境演變與生態建設實驗室，石家莊 050024)

0 引言

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

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

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

1 研究區概況與數據預處理

1.1 研究區概況

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

1.2 數據來源與處理

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

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

(1)

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

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

圖1 研究區概況圖、MODIS數據分幅及周邊氣象臺站

表1 MODIS數據產品

1.3 地面站點數據

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

表2 觀測點數據

圖2 氣象站點

1.4 其他數據

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

2 研究方法

2.1 隨機森林原理

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

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

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

LSTpre(EVI)=H(EVI)

(2)

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

圖3 降尺度示意圖

2.2 降尺度模型構建

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

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

圖4 技術流程圖

3 降尺度結果

3.1 降尺度模型篩選

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

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

3.2 降尺度模型RF的站點直接驗證

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

(3)

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

大興站、密云站、館陶站的降尺度LST與輻射四分量對比，平均誤差分別為3.86 K、2.37 K、1.96 K。地面氣象數據為點尺度的數據，對衛星數據進行驗證時理想的下墊面為勻質下墊面。其中，大興站為非勻質下墊面(旱地與果園)，因此在尺度匹配時出現系統誤差，導致降尺度地表溫度與地面氣象數據誤差較大；密云站、館陶站為勻質下墊面(旱地)，降尺度地表溫度與地面氣象數據誤差較小。

圖5 RF 模型的反演的LST與地面站點 LST散點圖

3.3 RF降尺度模型在不同下墊面的表現

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

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

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

表5 RF模型反演LST的評分

3.4 重建地表溫度缺失像元

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

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

4 結束語

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

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

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