李建明,馬燕飛,李仁杰,2,4,商國營,李明明
(1.河北師范大學 資源與環(huán)境科學學院,石家莊 050024;2.河北省環(huán)境變化遙感識別技術(shù)創(chuàng)新中心,石家莊 050024;3.邯鄲學院 地理系,河北 邯鄲 056005;4.河北省環(huán)境演變與生態(tài)建設(shè)實驗室,石家莊 050024)
陸地表面溫度(land surface temperature,LST)是陸地表面水熱平衡的重要分量,直接驅(qū)動地表與大氣之間能量交換。觀測與反演LST,對全球氣候變化評估、生態(tài)問題控制、水熱循環(huán)的進一步研究具有重要意義。同時,對農(nóng)業(yè)指導、旱情監(jiān)測、提高農(nóng)田水資源利用率做出了重要貢獻[1]。目前,LST已經(jīng)廣泛應(yīng)用到多種領(lǐng)域,如城市熱島[2]、土壤水估算[3]、蒸散發(fā)(evapotranspiration,ET)估計[4]等。但現(xiàn)有的地表溫度產(chǎn)品多為大網(wǎng)格尺度,難以滿足田塊尺度的地表熱通量研究[5]。MODIS熱紅外遙感地表溫度產(chǎn)品可以大面積獲取地表熱通量數(shù)據(jù),但是空間分辨率較低,在中緯度地區(qū)受云雨天氣影響,熱紅外地表溫度空間上不連續(xù)、缺失大量像元,實際應(yīng)用價值低[6-7]。近年來,星載熱紅外傳感器空間分辨率與時間分辨率不能共存的缺陷日趨明顯。在地表溫度反演領(lǐng)域引入LST降尺度方法,使獲取不同時空分辨率的地表溫度產(chǎn)品成為可能[8],因此,降尺度技術(shù)就成為連接大尺度LST與區(qū)域小尺度地表熱通量的紐帶。
地表溫度空間降尺度是在低空間分辨率像元值的基礎(chǔ)上,融合高空間分辨率的解釋向量,獲取更多的空間細節(jié)信息。降尺度根據(jù)算法原理的不同可以分為兩類。一類是基于物理機制的圖像融合降尺度方法,包括空間濾波分析法[9]、相關(guān)系數(shù)法(相關(guān)統(tǒng)計分析)[10]、主成分分析法等[11]。二類是基于尺度因子統(tǒng)計關(guān)系的降尺度方法[12]。Kustas等[13]建立陸地表面溫度與植被指數(shù)的統(tǒng)計回歸關(guān)系,融合高分辨率的航空數(shù)據(jù)(24 m)進行地球同步運行環(huán)境衛(wèi)星數(shù)據(jù)降尺度,得到田塊尺度的空間細節(jié)信息,并提出使用熱紅外遙感與NDVI進行TsHARP算法融合,對地表溫度反演;Ma等[14]使用NDVI在Kustas的TsHARP算法基礎(chǔ)上進行改進,對地表溫度進行降尺度。但是,傳統(tǒng)算法難以表征地表溫度與地表參數(shù)之間復雜統(tǒng)計關(guān)系。機器學習在地表溫度數(shù)據(jù)降尺度的探索、分析、預測和模型的建立與評估方面優(yōu)勢明顯,可以在眾多解釋向量中提取需要的信息,獲取高分辨率地表溫度影像,并重建熱紅外地表溫度的缺失像元。孟楊繁宇[15]使用機器學習的方法對黑河流域的蒸散發(fā)進行回歸、預測。受到啟發(fā),毛克彪等[16]嘗試將機器學習引入LST降尺度研究,使用神經(jīng)網(wǎng)絡(luò)算法對地表溫度進行反演,引入AMSR-E被動微波溫度數(shù)據(jù),并重建熱紅外數(shù)據(jù)在云雨天氣的數(shù)據(jù)缺失。Hutengs等[17]使用隨機森林算法對1 000 m分辨率的地表溫度進行降尺度,利用紅光、近紅外波段的地表反射率和數(shù)字高程模型DEM構(gòu)建最優(yōu)解釋向量集,建立簡單隨機森林降尺度模型(basic-RF),得到250 m分辨率的地表溫度。
本文以海河流域為研究區(qū),MODIS地表溫度產(chǎn)品為數(shù)據(jù)源,通過隨機森林算法對地表溫度進行降尺度,建立地表溫度與解釋變量(植被指數(shù)、空氣溫度、下行短波輻射)的降尺度模型。建立了精細圖像與粗圖像之間的非線性映射關(guān)系,選取植物生長月份進行地表溫度的回歸預測。探討隨機森林降尺度算法在不同下墊面的表現(xiàn),刻畫田塊尺度的地表溫度在海河流域的時空分布格局,為海河流域地表水熱循環(huán)研究提供數(shù)據(jù)參考。
海河流域位于112°E~120°E、35°N~43°N之間。東臨渤海灣,西倚太行山,南臨黃河,北接蒙古高原[18]。行政區(qū)劃包括京津冀大部、山西省東部、河南省北部及山東省東北部地區(qū)。流域總面積32.06×104km2,占全國總面積的3.3%。海河流域包括五大支流(潮白河、永定河、大清河、子牙河、南運河)和一個小支流(北運河)。海河流域地勢西高東低,以高原、山地以及平原等三種地貌為主。海河流域地處溫帶半濕潤半干旱大陸性季風氣候區(qū),流域年平均氣溫在1.5~14.0 ℃,年平均相對濕度在50%~70%之間,年平均降水量約為535 mm,年平均陸面蒸發(fā)量約為470 mm,水面蒸發(fā)量約為1 100 mm。近年來,海河流域氣溫呈明顯上升趨勢,水資源呈減少趨勢[19]。
本文主要使用MODIS產(chǎn)品中2~3級的V005版標準數(shù)據(jù)產(chǎn)品,包括地表溫度產(chǎn)品、植被指數(shù)NDVI產(chǎn)品,具體信息見表1。其中,地表溫度產(chǎn)品MOD11A1數(shù)據(jù)集的DN值(有效數(shù)據(jù))范圍為7 500~65 535,地表真實溫度值由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)方法進行投影的。如圖1所示,本文所需完整海河流域的MODIS數(shù)據(jù)是由四幅影像鑲嵌得到,需利用MODIS數(shù)據(jù)預處理工具MRT(MODIS reprojection tool)對這些產(chǎ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)品
“海河流域多尺度地表通量與氣象要素觀測數(shù)據(jù)集”[20]是國家青藏高原科學數(shù)據(jù)中心(https://data.tpdc.ac.cn/zh-hans/)發(fā)布的2008—2010年自動氣象站觀測數(shù)據(jù)。用戶可通過“寒區(qū)旱區(qū)科學數(shù)據(jù)中心”申請并下載該數(shù)據(jù)集。該數(shù)據(jù)集不僅可以定量揭示海河流域主要下墊面地表水熱通量交換特征,同時也為遙感估算地表蒸散量的驗證提供地面“相對真值”。如表2與圖2所示,數(shù)據(jù)來源主要包括三個自動氣象站點:海河流域密云站(山區(qū)林地)、館陶站(南部農(nóng)田)和大興站(中部城郊)。密云站下墊面是果園(李子樹和蘋果樹),館陶站下墊面是耕地(玉米、小麥、棉花),大興站下墊面是耕地(玉米、小麥),氣象站點下墊面性質(zhì)較為均一。自動氣象站的采集頻率為10 s,且10 min輸出一次。
表2 觀測點數(shù)據(jù)
圖2 氣象站點
本文還收集了用于地表溫度降尺度的海河流域部分地區(qū)的空氣溫度、下行短波輻射數(shù)據(jù)。用到的相關(guān)專題圖有:海河流域邊界矢量圖、海河流域土地利用(GlobeLand30 http://www.globallandcover.com)等。在具體應(yīng)用過程中還對這些專題圖件做了進一步的裁剪、重采樣及范圍匹配等處理。
隨機森林是用途廣泛的一種機器學習算法,隨機森林模型的本質(zhì)是對決策樹算法的改進。隨機森林模型是一個樹型分類器,輸入一個二維的矩陣,構(gòu)建沒有剪枝的回歸決策樹,每個葉子節(jié)點代表一個回歸結(jié)果,決定了單棵樹的生長過程,隨機森林的回歸模型采用單棵樹輸出結(jié)果的簡單平均得到。在估算數(shù)據(jù)的過程中,解釋變量的構(gòu)建對于目標數(shù)據(jù)相當重要。隨機森林是一個決策樹的集合,它的功能更多依賴于決策樹的功能[21],是在其基礎(chǔ)上的擴展,從訓練庫中提取實例進行判斷,根據(jù)屬性最終決定拿出哪個樣本來預測最優(yōu)結(jié)果。因此可以說,隨機森林是一個元估計器,適用于大量子樣本基礎(chǔ)上的決策樹分析,主要使用平均值來提高預測精度和控制過度擬合[22],最終的預測結(jié)果不會擺脫訓練樣本的范圍。本文構(gòu)建的模型本質(zhì)上屬于回歸問題,對于回歸問題隨機森林模型給出所有決策樹預測結(jié)果的平均。在地表溫度降尺度過程中,地表溫度與各地表參數(shù)(解釋向量)間呈非線性關(guān)系,而隨機森林降尺度模型對多元共線性不敏感。
本研究采用R語言中random forest數(shù)據(jù)包構(gòu)建隨機森林降尺度模型,訓練樣本為低空間分辨率(500 m、1 000 m)的遙感影像,地表溫度作為因變量,植被指數(shù)、空氣溫度、下行短波輻射作為解釋變量,以獲取更高空間分辨率的地表溫度。
如圖3所示,解釋變量集與地表溫度訓練出的某種映射關(guān)系構(gòu)成了隨機森林地表溫度降尺度模型,即在低空間分辨率影像(1 000 m)上建立地表溫度與地表參數(shù)(解釋變量)的統(tǒng)計關(guān)系,并假設(shè)這種統(tǒng)計關(guān)系不隨尺度的大小而變化。計算如式(2)所示。
LSTpre(EVI)=H(EVI)
(2)
式中:LSTpre表示地表溫度預測值;EV表示多個地表參數(shù)組成的最優(yōu)解釋變量集;I表示高空間分辨率影像(250 m);H表示地表溫度與最優(yōu)解釋變量集的回歸模型。
圖3 降尺度示意圖
隨機森林地表溫度降尺度由數(shù)據(jù)模型驅(qū)動,其本質(zhì)是以大量解釋向量樣本為基礎(chǔ),通過隨機森林模型方法訓練樣本,進而挖掘遙感圖像更加精細的空間細節(jié)。因此,解釋向量集的選擇十分重要。經(jīng)過實驗,本文將表3中的歸一化植被指數(shù)、空氣溫度與下行短波輻射構(gòu)成解釋向量集。本文地表溫度降尺度的核心是:隨機森林降尺度模型的建立與解釋變量集的選擇。降尺度模型構(gòu)建的具體步驟為:①MODIS LST數(shù)據(jù)與解釋變量數(shù)據(jù)的鑲嵌、重投影等預處理,并重分類為高分辨率(250 m)和低空間分辨率(500 m、1 000 m)的數(shù)據(jù)集;②將不同分辨率的解釋向量數(shù)據(jù)集、地表溫度數(shù)據(jù)集轉(zhuǎn)換成隨機森林算法可以識別的CSV文件;③將MODIS LST、解釋向量集作為降尺度模型的訓練數(shù)據(jù)進行訓練;④把高分辨率(250 m)的解釋向量數(shù)據(jù)輸入到訓練好的降尺度模型,獲取高分辨率(250 m)的地表溫度預測結(jié)果;⑤將降尺度模型預測結(jié)果CSV文件轉(zhuǎn)換回遙感圖像。技術(shù)線路如圖4所示。
表3 解釋變量集類型與名稱
圖4 技術(shù)流程圖
對隨機森林降尺度模型輸入?yún)?shù)進行調(diào)整,構(gòu)建500 m、1 000 m兩個不同分辨率層級的訓練樣本,使用均方根誤差(RMSE)、決定系數(shù)(R2)、偏差(bias)、相對精度(RA)對模型進行篩選,確定最優(yōu)模型算法。分別在500 m、1 000 m分辨率層級上使用隨機森林模型進行訓練。500 m與1 000 m分辨率層級模型的驗證結(jié)果如表4所示。在不同分辨率層級,相同的解釋向量組合的訓練集下,隨機森林降尺度模型的性能表現(xiàn)優(yōu)秀,可以使用更少的數(shù)據(jù)來描述解釋向量集與陸表溫度之間的非線性關(guān)系。在500 m和1 000 m隨機森林訓練模型中,隨機森林算法在1 000 m分辨率層級下表現(xiàn)出的空間異質(zhì)性更小,訓練模型表現(xiàn)更優(yōu),下一步采用1 000 m分辨率的訓練模型進行陸表溫度反演。
表4 隨機森林在500 m、1 000 m層級下模型反演LST的評分
地面觀測實驗的目的是為了獲取地面真實值,以便對遙感估算結(jié)果進行驗證,但是由于儀器觀測的不確定性和誤差的客觀存在,所以必須對觀測數(shù)據(jù)進行處理與質(zhì)量控制[23]。基于MODIS MOD11A1遙感數(shù)據(jù)來獲取地表溫度,降尺度后空間分辨率為250 m,這與自動氣象儀的觀測尺度是不匹配的。本文地表溫度觀測數(shù)據(jù)采用四分量輻射儀計算得到的結(jié)果。四分量輻射有效探測范圍分別是:館陶站15.6 m、大興站28 m、密云站30.76 m。一定程度上克服了站點尺度監(jiān)測范圍小的缺陷,計算如式(3)所示。
(3)
式中:σ為斯蒂芬-玻爾茲曼常數(shù)(5.67×10-8W/m2K4);ε為地表比輻射率;Rlu為大氣上行長波輻射;Rld為大氣下行長波輻射。
大興站、密云站、館陶站的降尺度LST與輻射四分量對比,平均誤差分別為3.86 K、2.37 K、1.96 K。地面氣象數(shù)據(jù)為點尺度的數(shù)據(jù),對衛(wèi)星數(shù)據(jù)進行驗證時理想的下墊面為勻質(zhì)下墊面。其中,大興站為非勻質(zhì)下墊面(旱地與果園),因此在尺度匹配時出現(xiàn)系統(tǒng)誤差,導致降尺度地表溫度與地面氣象數(shù)據(jù)誤差較大;密云站、館陶站為勻質(zhì)下墊面(旱地),降尺度地表溫度與地面氣象數(shù)據(jù)誤差較小。
圖5 RF 模型的反演的LST與地面站點 LST散點圖
為了定量描述隨機森林降尺度模型在不同下墊面的降尺度效果,如圖6(a)所示,本文在海河流域選取了五個不同下墊面性質(zhì)的子研究區(qū)。其中,子研究區(qū)A下墊面性質(zhì)為草地;子研究區(qū)B下墊面性質(zhì)為林地;子研究區(qū)C下墊面性質(zhì)為旱地;子研究區(qū)D下墊面性質(zhì)為旱地與建設(shè)用地;子研究區(qū)E下墊面性質(zhì)為林草交錯地。將五個子研究區(qū)降尺度結(jié)果升尺度至1 000 m,分別與MODIS原始地表溫度影像對比在不同下墊面性質(zhì)區(qū)域的降尺度效果并分析其空間分布特征。通過表5與圖6(c)可以看出,RF隨機森林方法在海河流域地表溫度降尺度上總體表現(xiàn)優(yōu)秀,這是由于RF隨機森林主要使用平均值來提高預測精度和控制過度擬合,可以更好地捕捉解釋向量與LST之間的函數(shù)關(guān)系,從而獲得高精度的降尺度結(jié)果。在植被覆蓋度高的草地、林地、林草交錯地帶降尺度效果R2、RMSE、bias優(yōu)于植被覆蓋度低的建設(shè)用地、旱地,降尺度的相對精度也更高,草地、林地、林草交錯地與原始溫度反演結(jié)果一致性更高。與植被覆蓋度高的子研究區(qū)相比旱地、建設(shè)用地RF隨機森林降尺度優(yōu)勢并不明顯。從影像的RMSE來看,建設(shè)用地達到3.894 K,說明解釋變量集(NDVI、空氣溫度、下行短波輻射)與LST的關(guān)系可能呈現(xiàn)出不穩(wěn)定的非線性關(guān)系。
從圖6(b)可以看出,在五個子研究區(qū)內(nèi),降尺度影像與MODIS原始影像的地表溫度在空間分布上一致,但總體來說降尺度影像的地表溫度要略低于MODIS原始影像。在降尺度地表溫度產(chǎn)品空間分布格局上,能夠合理地反映下墊面溫度變化,對地表溫度的空間分析具有重要意義。
注:A、B、C、D、E代表五個子研究區(qū);下標1表示降尺度地表溫度;下標2表示MODIS原始地表溫度。圖6 RF降尺度模型效果
表5 RF模型反演LST的評分
如圖7所示,通過RF|1 000 m|模型的降尺度的地表溫度數(shù)據(jù),融合250 m分辨率的無缺失像元解釋向量集得到的250 m海河流域地表溫度數(shù)據(jù),填補了MODIS MOD11A1地表溫度產(chǎn)品由云導致的地表溫度缺失像元。獲取海河流域地表溫度更精細的空間信息,得到精度更高的地表溫度產(chǎn)品,滿足海河流域地表溫度時空格局的研究。
注:①代表2008年第241天;②代表2010年第240天。圖7 海河流域降尺度結(jié)果
本文以海河流域為研究區(qū),探索了利用植被指數(shù)、空氣溫度、下行短波輻射等解釋向量對地表溫度空間降尺度的方法。通過對隨機森林算法在兩個空間尺度(500 m、1 000 m)下的訓練模型進行分析,隨機森林1 000 m算法表現(xiàn)優(yōu)于隨機森林500 m。將MODIS MOD11A1地表溫度產(chǎn)品與隨機森林模型降尺度的地表溫度在五個子研究區(qū)進行驗證分析,結(jié)果表明降尺度地表溫度的結(jié)果空間分布較為合理。將降尺度地表溫度與地面氣象站點的輻射四分量溫度進行直接驗證,在大興站、密云站、館陶站的平均誤差分別為:3.86 K、2.37 K、1.96 K。最終,構(gòu)建了隨機森林算法地表溫度降尺度模型。將MODIS MOD11A1地表溫度降尺度到250 m分辨率,并利用降尺度地表溫度填補MODIS地表溫度的缺失像元,避免了MODIS MOD11A1地表溫度產(chǎn)品空間分辨率較低、易受天氣影響、有效像元信息缺失嚴重、實際應(yīng)用價值低等問題。海河流域的降尺度地表溫度數(shù)據(jù)為數(shù)據(jù)稀疏區(qū)域的數(shù)據(jù)獲取提供科學參考。
本文總體上亮點有三點。一是將植被指數(shù)、空氣溫度、下行短波輻射組成的最優(yōu)解釋向量集引入到MODIS地表溫度降尺度,將地表溫度的空間分辨率由1 000 m降尺度到250 m,獲取到更加精細的地表溫度空間細節(jié)。二是在海河流域選出五個子研究區(qū),分析了在不同下墊面性質(zhì)下,隨機森林降尺度模型的效果。三是采用輻射四分量溫度對地表溫度降尺度結(jié)果進行直接驗證,一定程度上改善了驗證數(shù)據(jù)與降尺度地表溫度像元尺度不匹配的問題。
但文中仍存在一些不足需進一步研究,如通過遙感反演得到的地表溫度具有一定的方向(觀測角度)性,模型沒有考慮衛(wèi)星傳感器傾斜觀測角度造成的地表溫度差異。如何將傾斜觀測的地表溫度歸一到垂向觀測值,是未來研究需進一步考慮的內(nèi)容。