顏佳楠,陳 虹,姚光林,吳 驊
(1.中國科學院地理科學與資源研究所 資源與環境信息系統國家重點實驗室,北京 100101;2.中國自然資源航空物探遙感中心,北京 100083;3.山東省第八地質礦產勘查院,山東 日照 276826)
地表溫度(Land Surface Temperature,LST)是環境監測中的重要參數,可用于探測城市熱島效應[1],監測森林火災[2]和土地覆蓋變化[3],評估地表干旱程度[4]及管理水資源[5]。通過遙感衛星可以獲取到大范圍的LST,但是現有的衛星、遙感器受制于設備,無法獲取兼顧時間分辨率高和空間分辨率高的LST。例如中分辨率成像光譜儀(Moderate Resolution Imaging Spectroradiometer,MODIS)一天內可以獲取4次1 000 m分辨率的LST,而Landsat8衛星每16天獲取一次100 m分辨率的LST。LST遙感產品時空分辨率相互制約的局限性影響LST的應用。因此,許多學者嘗試針對MODIS數據,發展LST空間降尺度的方法,以期獲得時空分辨率較高的LST遙感產品。
目前LST空間降尺度方法主要有2種:一種是基于融合的方法;另一種是基于回歸的方法。基于融合的方法通過融合高時間分辨率和高空間分辨率的影像,預測高分辨率LST[6]。典型基于回歸的方法是Kustas等提出的DisTrad算法,其針對的是植被覆蓋區域,以NDVI為回歸核,獲取回歸核與低分辨率溫度之間的回歸關系,并將此函數關系應用于高分辨率回歸核與高分辨率LST[7]。Agam等人[8]基于TsHARP算法,將二次多項式回歸改為一次多項式回歸,把回歸核換作植被覆蓋度(FC),精度有所提高。然而,LST-FC(NDVI)的回歸關系在異質性地區表現較差[9]。針對這個問題,Essa等人[10]將回歸核改為SAVI,NDBI,NDBal,NDMI等15個光譜指標,將此方法推廣到異質性區域。Duan等人[11]以NDVI和DEM為回歸核,利用地理加權回歸進行LST降尺度。
為了提高模型的精度和適用范圍,有學者利用機器學習的方法獲取回歸核以及LST之間的關系,例如Gao[12]建立回歸樹獲取LST與短波光譜反射率之間的關系。Hutengs[13]提出了Basic-RF算法,用隨機森林的方法,獲取LST與數字高程數據、土地覆蓋類型、紅波段和近紅外波段的反射率之間的關系。Li[14]基于BasicRF提出了RRF算法,相比BasicRF,增加了回歸核種類,并對模型進行魯棒性驗證,在多個不同的研究區獲得了較高的精度。
考慮目前LST空間降尺度最優回歸核的組合形式以及其泛化能力尚不明確,本文嘗試選取波段反射率、遙感光譜指數、地形因子、地表覆蓋類型、經緯度以及大氣再分析數據6類回歸核,構建基于極端梯度提升樹(Extreme Gradrent Boost,XGBoost)算法的LST空間降尺度模型,同時開展對應的平行試驗,通過對比分析6類回歸核LST空間降尺度的精度,評估構建模型的泛化能力。
基于XGBoost算法的LST空間降尺度總體技術流程如圖1所示。

圖1 總體技術流程Fig.1 Flowchart of downscaling procedure
LST降尺度包括5個步驟:① 數據預處理。對圖像進行質量檢測、重投影、裁剪、重采樣以及空間聚合,獲取低分辨率(1 km)的回歸核數據和LST數據以及高分辨率(100 m)的回歸核數據和LST數據。② 溫度校正。由于待降尺度的MODIS的1 km分辨率LST產品是通過劈窗算法獲取的,而選作參考基準的Landsat-8衛星的100 m分辨率LST由單通道算法產生,搭載2個傳感器的衛星過境時間存在差異,為了更好地驗證LST空間降尺度的結果,本文將以Landsat-8衛星獲取的LST作為基準,對MODIS獲取的LST進行校正。③ XGBoost模型訓練。利用XGBoost獲取低分辨率回歸核與LST之間的回歸關系。④ 模型預測。用在低空間分辨率上訓練好的模型來預測高空間分辨率LST,加上低分辨率的殘差,得到溫度降尺度的結果。⑤ 降尺度結果的驗證。
XGBoost是由Chen[15]提出的一種提升樹(Tree Boosting)算法,該算法能夠應用于多種場景,在機器學習大賽和數據挖掘比賽中,多次被優勝隊伍選用[15]。Chen[15]提到XGBoost在梯度提升樹(Gradient Tree Boosting)的基礎上,進一步對模型做了優化。例如,在稀疏數據的處理上速度更快,并行和分布式計算等。這些優化使XGBoost能夠用最小的資源處理大量的數據。
在數據集D={(xi,yi)|i=1,2,…,n}(xi∈m,yi∈)中選取一個樣本(xi,yi),用K個可加性函數對樣本進行預測,結果如下:
(1)

XGBoost的目標函數由損失函數和正則化項組成,這2項分別表示模型的擬合效果和模型的復雜度。目標函數如下:
(2)

對損失函數進行二階泰勒展開,有:

基于回歸的LST降尺度方法的基本假設是在不同尺度下,LST回歸核的回歸關系不變[8],其關系式[16]為:

(3)


(4)

(5)
本文通過分析比較前人研究,共篩選出6類不同的回歸核,具體情況如表1所示。地表反射率表征地表反射的太陽短波輻射能量大小,與地面狀態、太陽高度角相關。光譜指數、地表類型數據提供地表狀態信息,包括植被覆蓋度、干旱狀態及建筑信息等,提高模型在異質性區域的適用性。地形數據反映不同地形地貌對輻射的影響,提高模型在地形起伏區域的適用性。大氣再分析數據反映大氣生物物理量,淺層土壤含水量受淺層土壤溫度的影響,且淺層土壤溫度、近地表氣溫與地表溫度有相關性。空氣流動的速率與植物蒸騰作用相關,間接對地表溫度產生影響。經緯度提供空間變化信息。

表1 回歸核的信息
為了比較不同回歸核情況下LST降尺度精度的差異,對6類回歸核進行不同組合,共得到12組,具體分組情況如表2所示。

表2 分組情況
依照地理位置分布范圍廣、地表類型多樣以及Landsat8圖像晴空無云且質量較高3個標準,以2019年為待選年份,在中國范圍內選出了7景滿足條件的Landsat8數據,并將其作為降尺度的研究區(83°45′E~116°17′E,33°N~44°30′N),研究區的時空信息如表3所示。

表3 研究區時空信息
本研究所用的Landsat8 C2L2數據來源于美國地質調查局網站(https:∥earthexplorer.usgs.gov/)。Landsat8數據集包含30 m分辨率的可見光波段、近紅外波段以及100 m分辨率的熱紅外波段。分別對原始Landsat8可見光-近紅外反射率產品以及熱紅外LST產品進行裁剪、重采樣(雙線性插值法)、空間平均聚合,得到100 m分辨率和1 km分辨率的LST以及相應的降尺度回歸核。
低分辨率的LST選擇分辨率為1 km的MOD11A1數據[17](https:∥lpdaac.usgs.gov/products/mod11a1v006/)。對MOD11A1數據進行重投影、裁剪,并利用質量文件刪除質量較差以及有云的點。
數字高程數據來源于航天飛機雷達地形測繪使命(Shuttle Radar Topography Mission,SRTM)[18](https:∥srtm.csi.cgiar.org/srtmdata/),分辨率為90 m,經過重投影、重采樣(雙線性插值法)、裁剪及空間平均聚合后,作為100 m分辨率和1 km分辨率的回歸核。
地表覆蓋類型數據來源于GlobeLand30數據集[19](https:∥www.webmap.cn/commres.do?method=globeIndex),分辨率為30 m。同樣對其進行重投影、重采樣(雙線性插值法)、裁剪、空間平均聚合的數據預處理工作,得到100 m分辨率和1 km分辨率的數據。
大氣再分析數據取自歐洲中期天氣預報中心(ECMWF)全球氣候大氣再分析資料的第五代產品——ERA5[20](https:∥cds.climate.copernicus.eu/),其空間分辨率為0.25°。從中選取淺層土壤溫度、淺層土壤體積含水量、離地2 m附近大氣溫度、離地10 m附近徑向風速、離地10 m附近橫向風速,5個數據作為代表。首先獲取與衛星成像時間最近2個時刻的大氣再分析數據,對大氣再分析數據進行重投影、裁剪以及雙線性插值處理,得到100 m分辨率和1 km分辨率的數據;再對相鄰2個時刻的數據進行線性插值,得到與衛星成像時間相同的大氣再分析數據。
此外,由于Landsat8獲取的LST與MOD11A1獲取的LST之間存在算法和觀測時間之間的差異,對MOD11A1的溫度產品進行了溫度校正。主要操作是先對Landsat8獲取的LST進行空間平均聚合操作,得到1 km分辨率的LST(Landsat8),再以MOD11A1-LST為自變量,Landsat8-LST(1 km)為因變量做回歸,對原始待降尺度的MOD11A1-LST進行了溫度的校正。7個研究區溫度校正過程的散點圖如圖2所示(圖2(a)~圖2(g)分別是7個研究區的校正結果),LST校正的均方根誤差(RMSE)在0.6~1 K,表明溫度校正后的結果可用于LST降尺度結果的檢驗。

圖2 MOD11A1溫度校正的散點圖



表4 12組平行試驗LST降尺度算法性能統計表
為了更好地展示各研究區降尺度的性能,本文挑選了表征降尺度精度的RMSE以及表征VIF作為評價指標,統計結果如表5所示,其中研究區1~7用A~G標識,分組用G1~G12表示,RstdV和VstdV分別表示7個研究區RMSE和VIF統計指標的標準差。

表5 各研究區的均方根誤差和視覺信息保真度
由表4和表5可以看出,group1(回歸核為地表反射率)的RMSE結果較差,表現為降尺度的總體精度較差,且各個研究區的離散度較大,標準差為0.91;剩余group的降尺度總體精度差別不大,在2 K左右,但group3(回歸核為地表反射率、光譜指數)在各個研究區的表現略微差,體現在精度存在部分差異,標準差在0.7左右。對于視覺信息保真度VIF,group1(回歸核為地表反射率)、group2、group3、group4(回歸核為光譜指數、地形)、group6(回歸核為光譜指數、地表類型)、group7、group12(回歸核為地表反射率、經緯度)能夠更好地保持圖像的紋理。
為了進一步目視展示降尺度的結果,本文選取了研究區1、研究區2和研究區7,繪制了不同分組回歸核降尺度前后的對比圖,如圖3所示。
由圖3可以看出,以降尺度的清晰度和LST的空間分布作為評價指標,結合表4和表5統計結果,group2、group4、group7相對更優,LST降尺度RMSE可控制在2 K,LST回歸殘差控制在0.5 K,而視覺信息保真度VIF可超過0.07。由于這3個分組中都包含了光譜指數回歸核,因此光譜指數對于LST的降尺度更加重要。在構建回歸模型時,引入光譜指數回歸核,能夠有效控制LST降尺度的精度。
為探究利用單個研究區進行訓練的模型推廣至其他研究區后降尺度精度、保真度指標以及回歸模型適用性的變化,即探尋降尺度模型的泛化能力,選擇3組(group2,group4,group7),針對研究區1,2和7,分別統計了選定研究區不同分組回歸核條件下的降尺度的性能指標,如表6所示。

表6 不同組別的回歸核在不同研究區交叉驗證的統計結果

圖3 選定研究區降尺度前后的對比圖Fig.3 Results of comparison for studied areas before and after LST downscaling
由表6可以看出,對于group2、group4和group7分組回歸核,某一研究區數據訓練得到的XGBoost模型應用于本研究區或者其他研究區時,LST降尺度的RMSE變化不大,二者差值小于0.5 K。對于視覺信息保真度VIF,某一研究區數據訓練得到的模型應用于其他研究區時,VIF將降低,會造成LST降尺度結果的失真,如模糊或者細節丟失。對于LST回歸模型的殘差,殘差變化較大,本地訓練模型應用于本地時,殘差值較小,然而用于其他區域的時候,殘差能增大到10 K以上。雖然LST降尺度的RMSE能保持在同一水平,但殘差值仍然較大,說明通過局部數據訓練的XGBoost模型,回歸泛化能力不夠,某個研究區訓練得到的模型,應用于更大范圍的研究區時可能會帶來較大的誤差。為了提高XGBoost模型的回歸泛化能力,需要保證訓練數據的代表性。
本文選取7個研究區,并對地表反射率、光譜指數、地形因子、大氣再分析數據、經緯度、地表類型6種回歸核進行組合實驗。利用基于XGBoost算法的LST降尺度模型,在不同回歸核組合的情況下,將1 km分辨率的MODISLST空間降尺度成100 m分辨率的LST,并將估計值與Landsat-8反演的100 m分辨率的LST進行了分析比較,討論了不同回歸核組合條件下LST降尺度的性能差異。本文得到的主要結論如下:
① 通過12組回歸核的實驗比對發現, group2,group4,group7作為回歸核時,可以得到精度低且質量好的預測圖像,3組的RMSE在2 K左右,LST降尺度殘差控制在0.5 K左右,VIF可超過0.07。
② 由于選定的最優回歸核組合(group2,group4,group7)中均包含了光譜指數,因此光譜指數是LST降尺度的關鍵因素。在實際應用中,可根據其他回歸核(如地形、經緯度)的獲取情況選擇具體的回歸核組合。從當前研究結果來看,引入更多的回歸核會部分提高LST降尺度的精度,但在一定程度上增大了降尺度實施的復雜度。
③ 基于XGBoost算法的LST降尺度模型的泛化能力仍舊不夠,雖然LST降尺度的RMSE可以得到一定的保證,但這是通過殘差糾正來實現的,實際上局部數據訓練的XGBoost回歸模型自身無法有效從某一區域推廣到另一區域。為了構建適用大范圍的LST降尺度模型,確定XGBoost回歸模型需要選擇更具代表性的訓練數據。
在本文研究的基礎上,后續可嘗試對降尺度模型進行改進,減少降尺度殘差對結果的影響。同時將考慮組內回歸核的相關性,對組內回歸核數據進行主成分分析,并對回歸結果進行比較分析。此外,也將重點考慮回歸窗口的影響,采用全局和局部等不同策略來完善降尺度模型,以期得到一個泛化能力更強的LST降尺度模型。