喀斯特地區地表溫度空間降尺度方法初探

尹 枷 愿，蔡 宏*，田 鵬 舉，唐 敏

(1.貴州大學礦業學院，貴州 貴陽 550025；2.貴州省生態氣象和衛星遙感中心，貴州 貴陽 550025)

0 引言

地表溫度(LST)是大氣與地表能量循環中最重要的參數之一[1]，被廣泛應用于水文平衡評估、全球變暖研究、城市熱島效應評估和地表蒸散量計算等領域[2-7],然而，受熱紅外傳感器成像條件的影響，地表溫度數據產品常因時空分辨率相互制約而無法滿足實際應用需求。為此，不少學者提出地表溫度降尺度方法以解決該問題，其實質是將其他高分辨率數據信息整合到熱紅外波段上，以此提高地表溫度分辨率。常用的地表溫度降尺度方法主要分為基于圖像融合降尺度方法和基于尺度因子的統計回歸降尺度方法[8]，后者著重研究影響LST變化的物理因素，因操作簡單、精度高等而得到廣泛應用。

早期統計回歸降尺度研究主要探究LST與歸一化植被指數(NDVI)[9]、植被覆蓋度(VC)[10-12]等植被指數之間的線性關系，并成功應用于地形均勻或同質地區，但在復雜地區的效率不高[13]。因此，不少學者加入更多凸顯地形變化和地表類型的尺度因子，例如：Maeda利用海拔、NDVI與LST之間的線性回歸關系進行山區降尺度研究[14]；Duan等通過地理加權回歸方法建立NDVI和DEM之間的非平穩關系，并成功將河北省華北平原與內蒙古高原過渡地帶的LST從1 km降尺度至90 m[15]。盡管傳統的統計回歸方法取得了一定成效，但考慮到尺度因子與LST之間復雜的非線性關系，一些學者建議將貝葉斯[16]、人工神經網絡[17]、支持向量機[18]和隨機森林(RF)[19]等機器學習模型用于地表溫度降尺度。其中，RF具有實現簡單、不易過擬合的特點，能更好地模擬LST與地形、地貌間的關系。例如：Hutengs等利用DEM、土地覆蓋圖、紅光與近紅外波段的反射率同LST建立RF模型，在地中海東部高程變化大的植被地區，將LST從1 km降尺度到250 m[19]；Wu等利用地表反射率、光譜指數、地形因子和土地利用圖，對以自然地表為主的西班牙塞哥維亞佩尼亞羅拉山區和北京平原城區進行了降尺度研究[20]。然而，目前考慮喀斯特山區地形地貌異質性特點的高精度地表溫度降尺度鮮有研究。因此，本文以遵義市西北部和貴陽市以南地區為研究區，選擇適合喀斯特山區的尺度因子，以MODIS 第31、32波段的輻射亮度為因變量，建立一種適合喀斯特地區的隨機森林(Karst Random Forest,KRF)模型，并以Landsat-8反演的地表溫度數據為驗證數據進行精度評定。

1 研究區及數據

考慮到喀斯特地區地形起伏度相差較大，本研究選擇兩個有代表性的區域(圖1)，以充分探究KRF模型的適用性。區域a位于遵義市西北部，為喀斯特山區向平原過渡的盆地邊緣地帶，處于副熱帶東亞大陸季風區內，屬于中國亞熱帶高原季風濕潤氣候，區域內地勢起伏較大(最高海拔為1 730 m，最低海拔為316 m，高差為1 414 m)，地表覆蓋類型主要為自然地表及小部分建設用地和水體；區域b位于貴陽市以南，為典型喀斯特山區城市，最高海拔約為1 660 m，最低海拔約為991 m，高差為669 m，氣候類型與區域a相同，地表覆蓋類型主要有建設用地、水體、裸土、裸巖等。

圖1 研究區示意

研究數據(表1)包括：1)MODIS數據，來源于NASA(https://ladsweb.modaps.eosdis.nasa.gov/)，選用其MOD021KM和MOD11A1 產品，軌道號為h27v06，成像時間分別為2016年7月26日和2018年3月4日，利用NASA發布的HEGTool軟件將該數據重投影為WGS-84/UTM。2)Landsat-8 OLI/TIRS數據，來源于地理空間數據云(http://www.gscloud.cn/)，區域a條帶編號為128，行編號為40，選取與MODIS同一天(2016年7月26日)的數據；區域b條帶編號為127，行編號為42，由于貴陽地區云霧較大，無法獲取同一天的數據，通過對歷史氣象站數據查證，選取成像時間與MODIS日期相近的2018年3月3日的Landsat-8影像作為參考影像。3)高程數據，來源于地理空間數據云的ASTER GDEM V2，時間為2015年1月6日，條帶編號為105、106，行編號為26、28。以上影像利用ENVI進行輻射定標、大氣校正、鑲嵌和裁剪等預處理。

2 喀斯特地區地表溫度空間降尺度方法

喀斯特地區地表溫度空間降尺度流程包括尺度因子選取、RF熱紅外降尺度模型構建、熱紅外地表溫度反演和精度評價四部分(圖2)。

2.1 尺度因子選取

尺度因子會影響降尺度方法的精度，本研究針對喀斯特地區的特點，添加了常規尺度因子和喀斯特地表特征尺度因子。

表1 研究數據信息

圖2 喀斯特地區地表溫度降尺度流程

2.1.1 常規尺度因子 著重考慮表征土地覆蓋的相關指數，篩選后的參數包括Landsat-8 OLI的紅、藍、綠、近紅外和短波紅外波段的反射率以及Landsat-8 OLI反演的遙感指數：表征植被的比值植被指數(RVI)、土壤調整植被指數(SAVI)、歸一化植被指數(NDVI)、植被覆蓋度(VC)，表征水體的改進歸一化水體指數(MNDWI)、歸一化干旱指數(NDDI)，表征城鎮的城市指數(UI)、建筑用地指數(IBI)。

2.1.2 喀斯特地表特征尺度因子 Deng等在研究喀斯特地區土地利用、NDVI與LST的關系時，發現低植被覆蓋的裸巖與裸土區域存在異常高溫[21]，孟鵬燕等發現地形是影響LST的重要因素[22]，因此，除常規尺度因子外，本文在地貌上選擇裸土指數(BSI)[23](式(1))和裸巖率(BRp)[24](式(2))，地形上選擇由ASTER GDEM數據得到的高程因子，包括高程、坡度、坡向、山體陰影。所有尺度因子的空間分辨率均為30 m，通過平均聚合到1 km與MODIS數據進行降尺度建模。

(1)

BRp=1-BSp-VC

(2)

式中：ρmir為中紅外波段反射率，本文用Landsat-8的短波紅外波段1代替;ρred、ρnir、ρblue分別為紅光、近紅外、藍光波段的反射率;BSp為土壤裸露率；VC為植被覆蓋度。

2.2 基于RF算法的熱紅外降尺度模型構建

貴州喀斯特地區地表形態特殊，景觀破碎且異質性強，復雜的地表環境使尺度因子對地表溫度的響應關系復雜，因此，本研究考慮尺度因子與熱紅外輻射亮度之間的非線性關系，引入隨機森林(RF)算法進行建模。該算法由大量不相關的決策樹構成，通過bootstrap重抽樣方法，首先將原始訓練樣本中抽取的多個樣本自主合并，生成訓練樣本合集，然后根據自助樣本集生成多個決策樹并組成隨機森林，最終結果由所有決策樹平均而得，具有精度高、對多元線性不敏感的特點。

基于尺度因子的LST降尺度方法的本質是：不同尺度下LST和尺度因子之間的定量關系保持不變，即在低空間分辨率尺度下建立的關系模型，在高空間分辨率尺度下依然適用。降尺度模型的表達式為：

(Bi)C=f((Ak)C)

(3)

ΔC=(Bi)O-(Bi)C

(4)

(Bi)F=f((Ak)F)+ΔC

(5)

式中：i=31、32，分別代表MODIS第31、32波段；(Bi)C為模型計算得到的第31、32波段低分辨率輻射亮度；(Ak)C為通過平均聚合得到的低空間分辨率的尺度因子，k為尺度因子數量；f(·)為RF回歸映射關系；ΔC為原始影像與預測影像之間的殘差；(Bi)O為第31、32波段原始輻射亮度；(Bi)F為模型計算的第31、32波段高分辨率輻射亮度；(Ak)F為高分辨率的尺度因子。

本文中RF模型在Python平臺上實現，其中關鍵參數為決策樹的數量(n_estimators)和樹的最大特征數量(max_features)，通過交叉驗證、遍歷優化，最終確定區域a第31波段RF模型的n_estimators=290，max_features=17，第32波段RF模型的n_estimators=375，max_features=13；區域b第31波段RF模型的n_estimators=390，max_features=11，第32波段RF模型的n_estimators=391，max_features=14。

2.3 熱紅外地表溫度反演

現有LST反演算法主要包括單通道算法、劈窗算法和多通道多角度算法3類。針對不同的傳感器，選擇不同的反演方法能有效提高反演精度。

2.3.1 MODIS地表溫度反演 MODIS第31、32波段的中心波長接近于AVHRR的第4、5通道，適宜用劈窗算法進行LST反演，因此本文選用覃志豪等的劈窗算法[25]，對降尺度后的第31、32波段輻射亮度進行反演，得到LST。反演公式為：

Ts=A0+A1T31-A2T32

(6)

(7)

(8)

(9)

Ci=εiτi(θ)

(10)

Di=[1-τi(θ)][1+(1-εi)τi(θ)]

(11)

式中：T31、T32分別為MODIS第31、32波段的亮溫，通過普朗克公式計算；Ai、Ci、Di為參數；a31、b31、a32和b32根據波段特征確定，在地表溫度0～50 ℃范圍內，分別取a31=-64.60363，b31=0.440817，a32=-68.72575，b32=0.473453；εi為地表比輻射率，可由土地利用圖和NASA提供的比輻射率計算；τi(θ)為大氣透過率，可通過MOD021KM數據計算得到。

2.3.2 Landsat-8地表溫度反演 Landsat-8 TIRS的第11波段不穩定，不宜反演地表溫度，因此本文使用改進的單窗算法[26]對第10波段進行地表溫度反演(式(12)-式(14))，用于降尺度驗證。由于Landsat-8熱紅外波段的原始分辨率為100 m，所下載的為NASA重采樣到30 m的數據，并與30 m的OLI數據進行了配準[27]，因此，本文參照張文奇等的方法，利用聚合平均法對TIR數據重采樣后再驗證[28]，即將Landsat-8第10波段從30 m重采樣到100 m，再進行溫度反演。

Ts=[a10(1-C10-D10)+(b10(1-C10-D10)+C10+D10)T10-D10Ta]/C10

(12)

C10=τ10ε10

(13)

D10=(1-τ10)[1+(1-ε10)τ10]

(14)

式中：a10和b10是第10波段的普朗克回歸系數，分別取-62.7182和0.4339。

2.4 精度評價指標

為評估本文適合喀斯特地區的隨機森林(KRF)方法的精度，分別與目前最常用的TsHARP算法(基于NDVI與LST間線性關系的經典降尺度算法)[9]和最新的MTVRF模型(基于反射率、非喀斯特光譜指數、土地覆蓋圖及高程等18個因子與LST間非線性關系的RF降尺度算法)[20]進行對比，并分別從總體和各地類兩方面對降尺度結果進行定性和定量分析。定性方面主要對比3種降尺度方法得到的LST和Landsat-8反演的LST在空間分布特征上的差異；定量方面主要通過均方根誤差(RMSE)、平均絕對誤差(MAE)衡量降尺度結果與真實值間的偏差，其值越小，說明預測值與實際值的差值越小，降尺度精度越高；應用線性回歸決定系數(R2)評定回歸模型的擬合優度，R2越接近1，說明回歸模型的擬合程度越高。具體計算公式如下：

(15)

(16)

(17)

3 結果與討論

3.1 熱輻射亮度和尺度因子的重要性評價

為探究熱輻射亮度和尺度因子對降尺度的影響，首先將選定的20個尺度因子分別與MOD11A1的第31、32波段的熱紅外輻射亮度進行RF建模，交叉驗證5次，得到區域a 尺度因子與MOD11A1數據的R2為0.477，小于第31波段R2(0.552)和第32波段R2(0.513)，區域b 尺度因子與MOD11A1數據的R2為0.556，小于第31波段R2(0.642)和第32波段R2(0.634)，這與Agathangelidis等[29]利用輻射亮度與尺度因子建模后精度更高的結論一致。因此，直接利用熱紅外輻射亮度進行降尺度能夠有效減少誤差，提高精度。

為進一步分析尺度因子對降尺度的影響，通過尺度因子與輻射亮度的RF模型，得到各尺度因子對不同波段輻射亮度的貢獻度(圖3)。由圖3可知，DEM在兩區域對不同波段的貢獻度均較高，說明山區地表溫度空間分布受地形變化引起的熱環境差異影響，地形效應顯著[30]；此外，區域a的BSI、RVI對不同波段輻射亮度的貢獻度均較高，這可能是因為該區域為地形高差大的自然地表，地表溫度空間分布主要受地形因素和表征自然地表的遙感參數控制；區域b中SWIR2的貢獻度較高，這可能是因為該區域為城區，總體地表溫度較高，而SWIR2波段對高溫目標較敏感[31]。

圖3 尺度因子貢獻度評估

整體而言，不同尺度因子對熱紅外輻射亮度均有貢獻，表明本文所選的尺度因子適用于研究區，但不同區域因地形地貌不同，因子的貢獻度有所不同，而RF模型會根據所添加的因子與熱紅外輻射亮度的關系自行調節權重，因此，加入更多適合喀斯特地區的參數進行降尺度后，能更好地模擬出喀斯特地區不同地形地貌對熱紅外輻射亮度的影響，使降尺度模型在不同區域均表現出較高的精度，從而提高模型的普適性。

3.2 不同喀斯特區域降尺度精度評價

利用3種方法在兩個代表性喀斯特區域分別對MODIS數據進行降尺度，結果如圖4、圖5(彩圖見附錄1)所示，并以Landsat-8反演的100 m LST為真實值進行對比，發現3種方法均可將LST的空間分辨率從1 km提高到100 m。其中，TsHARP算法只考慮植被與LST的線性關系，在植被與水體和城鎮的邊界處噪聲較多，效果最差；MTVRF模型和本文KRF方法為多尺度因子降尺度模型，在兩區域精度都較高，但KRF方法效果更優，尤其是在中溫區和高溫區的空間分布上可顯示更多細節，更加精細地描繪出喀斯特地區熱特征空間分布的差異性。

圖4 區域a總體降尺度結果比較

圖5 區域b總體降尺度結果比較

為更客觀地比較3種降尺度方法在不同喀斯特地區的表現，分別計算各方法的MAE、RMSE和R2(表2)。可以看出，對于區域a，KRF方法的性能較好，其MAE和RMSE分別為2.1146 K、2.4652 K，較TsHARP算法分別減少了0.4439 K和0.6204 K，較MTVRF模型分別減少了0.1419 K和0.1302 K，R2較TsHARP算法有一定的提高，但較MTVRF模型有所下降。對于區域b，KRF方法的MAE為1.1363 K、RMSE為1.4521 K，較TsHARP算法分別減少了0.5084 K 和0.6953 K，較MTVRF模型分別減少了0.214 K 和0.2928 K，R2為0.7103，較TsHARP算法和MTVRF模型分別提高了0.3805和0.1556。綜上可知，本文KRF方法在喀斯特地區效果最好，說明引入與喀斯特相關的指數有助于降尺度，且利用MODIS第31、32波段進行建模能夠最大限度避免異常值。

表2 喀斯特地區不同地形上的降尺度結果精度

另外，自然地表區域a的精度小于城市區域b的精度，這與Wu等[20]的研究結果相同，可能是因為地形變化對LST的影響更大。因此，利用KRF方法雖能明顯提高喀斯特地區的地表溫度降尺度精度，但在地形起伏度相對較小的喀斯特山區城市反演精度更高，在地形起伏度較大區域精度略低。

3.3 喀斯特地區各地類降尺度精度評價

為了更清楚地顯示喀斯特地區局部降尺度成果，分別從區域a和區域b 4種地類(建設用地、植被、裸土、水體)中截取最具代表性的局部進行對比(圖6-圖9，彩圖見附錄2)，可以看出3種降尺度方法在不同地類區域的準確性不同，但總體上KRF方法在4種地類上的表現均優于其他算法。TsHARP方法僅在同質的植被區表現較好[18]，在其余地類的LST結果中均出現低估現象；MTVRF方法引入了較多的尺度因子，雖然在水體、植被區域表現較好，但在建設用地和裸土上表現較差，與Wu等[20]的研究結果不同，這可能是因為喀斯特地區建設用地和裸土分布較為破碎，導致這兩個區域與其他地物的交界處出現異常值。

進一步統計4種地類的評價指標(表3)，可以看出：KRF方法的MAE比MTVRF模型和TsHARP算法分別降低了0.2069～0.4907 K和0.3596～3.6282 K，RMSE比MTVRF模型和TsHARP算法分別降低了0.2818～0.7051 K和0.4129～4.184 K，R2比MTVRF模型和TsHARP算法分別提高了0.1014～0.1965和0.1329～0.4688，說明本文方法在4種地類中的精度更高；另外，3種方法在植被和水體區域效果較好，但在裸土和城鎮區域出現較大偏差，這可能是因為裸土和城鎮區域異質性較強，說明精度與地表異質性相關[32]。

圖6 喀斯特地區水體降尺度結果比較

圖7 喀斯特地區植被降尺度結果比較

圖8 喀斯特地區建設用地降尺度結果比較

圖9 喀斯特地區裸土降尺度結果比較

表3 3種方法在4種地類的降尺度結果精度

4 結論

本文從熱紅外輻射亮度出發，針對喀斯特地區的特征，添加合適的尺度因子，利用RF模型將空間分辨率為1 km的MODIS 第31、32波段輻射亮度降尺度至100 m，通過劈窗算法對100 m的熱紅外輻射亮度進行反演，最終構建了適合喀斯特地區的隨機森林(KRF)模型，實現了喀斯特地區的地表溫度空間降尺度。從喀斯特地區起伏度不同的兩個區域和4種不同地類兩種尺度上，將本文方法與MTVRF模型、TsHARP算法進行對比，得到如下結論：1)在不同起伏度的兩個喀斯特地區，本文方法均能較好地模擬LST的熱特征，且更適用于喀斯特地形變化較為平緩的山區城市。2)在不同地類降尺度方面，本文方法表現最好，特別是添加了BSI、BRp和高程因子與熱輻射亮度的降尺度研究，拓展了傳統降尺度模型的應用范圍，顯著提高了喀斯特地區各地類的降尺度效果。整體而言，本文方法較MTVRF模型和TsHARP算法在定性和定量上都得到明顯改善，能有效提高喀斯特地區地表溫度空間分辨率，但針對地形變化更大、異質性較強的區域，降尺度結果還有改進的空間。