基于隨機森林算法對青藏高原TRMM降水數據進行空間統計降尺度研究

徐彬仁， 魏瑗瑗

(1.中國科學院遙感與數字地球研究所，北京 100101； 2.中國科學院大學，北京 100049)

0 引言

水循環促進了自然界的物質運動和能量交換，對氣候的形成與變化產生了深刻影響。流域降水量是影響流域水循環最重要的因素。有“亞洲水塔”之稱的青藏高原是亞洲大江大河的發源地，印度河、恒河、雅魯藏布江、長江和黃河均源自青藏高原[1-2]。因此，青藏高原的流域的降水資料對我國乃至全球的水文、氣象、生態和農業等領域的發展具有重要的研究意義。與氣象站點降水資料相比，衛星資料空間覆蓋連續，有助于解決研究區站點數量不足、分布不均勻的問題。但是，由于探測平臺的限制，衛星遙感降水產品的分辨率仍然無法滿足流域尺度研究的要求。為進一步了解青藏高原流域尺度上降水的時空分布特征，需要開展青藏高原遙感降水數據的空間降尺度方法的研究。

目前，空間降尺度方法主要有動力降尺度和統計降尺度2種。空間動力降尺度是借助全球環流模式(global climate models, GCMs)和嵌套區域氣候模式(regional climate models, RCMs)來提高氣象要素的空間分辨率。該方法不受觀測資料的影響，但其計算量大、模型不易構造，獲得的氣象要素的空間分辨率不能滿足區域尺度的要求[4]。統計降尺度充分考慮到局地氣候不僅以大尺度氣候為背景，且顧及下墊面特征的影響，利用下墊面特征信息建立大尺度和小尺度氣候變量的聯系[5]。與空間動力降尺度方法相比，統計降尺度方法構造的模型更加靈活多樣，引入較高空間分辨率的區域下墊面特征變量，能大大提高遙感降水資料的空間分辨率[6]。2009年Immerzeel團隊以歐洲南部伊比利亞半島為例，基于降水與植被的關系，采用熱帶降水測量計劃衛星(Tropical Rainfall Measuring Mission, TRMM)降水數據與SPOT -VEGETATION 歸一化植被指數(normalized difference vegetation index, NDVI)數據，建立了的指數回歸降尺度模型，最終將TRMM降水數據降尺度為1 km分辨率。基于Immerzeel的研究，2011年賈紹鳳研究團隊開展了柴達木盆地降水降尺度研究，同時引入了植被數據SPOT-VEGETATION NDVI和航天飛機雷達地形測繪任務地形數據(Shuttle Radar Topography Mission，SRTM) 數字高程模型(digital elevation model，DEM)因子作為流域降水的影響因子，建立了多元線性回歸模型，也得到了1 km空間分辨率的降水數據[8]。與Immerzeel團隊的研究相比，后者引入DEM變量代表地形因子參與建模，并對降水的降尺度結果進行了誤差分析，進一步提高了降水資料的準確性。此外，非線性時間序列可以提高降水的預測精度[16]。在諸多非參數統計回歸模型中，隨機森林算法在分類和預測方面對自變量的多元共線性不敏感，可以同時輸入多個影響因子，在很大程度上解決了過度擬合的問題等。因此，本文基于Immerzeel和賈紹鳳研究團隊的研究和隨機森林算法的優勢，選擇NDVI,DEM,坡度、坡向和經緯度信息，針對青藏高原流域的長時間序列TRMM遙感降水數據，開展了空間統計降尺度分析研究。

1 研究區概況

青藏高原平均海拔4 500 m，整體面積達250×104km2，是世界上平均海拔最高的地區。獨特的地理環境使其對亞洲地區以及世界氣候具有重要的影響[1]。該地區地面氣象觀測站點主要集中于東部和南部，極少數分布在高原西部和北部，存在氣象資料空間分布不勻及稀缺現象。已有站點大多安置于低海拔的交通可抵達區域，并且觀測站測量的氣象要素受其周圍區域小氣候影響較大[14]。為彌補青藏高原地區傳統觀測數據稀缺、受局地氣候影響和空間分布不連續等不足，本文采用TRMM3B43降水數據開展了青藏高原地區降水量的空間統計降尺度研究。

受孟加拉灣暖濕空氣的影響，雅魯藏布江大拐彎區域年降水十分充沛； 高原的腹地湖泊數量較多，湖泊集聚區的降水分布普遍高于周圍地區； 而受喜馬拉雅山體影響，來自印度洋的水汽輸送受到阻擋，喜馬拉雅山北坡的降水明顯少于周圍地區。總體而言，青藏高原降水空間分布呈現自東南向西北遞減、自南向北逐漸減少的趨勢[14]。

2 數據源及其處理

2.1 TRMM衛星降水數據與氣象站降水數據

TRMM是世界上第一顆搭載測雨雷達的衛星，攜帶了微波成像儀、可見光和紅外掃描儀[15]。TRMM3B43產品綜合了4類相互獨立的降水數據，包括微波、近紅外等傳感器融合估算數據，以及美國國家海洋、大氣管理局和全球降水氣候中心的降水雨量計分析數據等[14]，是衛星降水數據與其他降水數據聯合反演的最佳降水產品。本文選用2000—2012年間空間分辨率為0.25°×0.25°的TRMM3B43日降水產品，對青藏高原地區進行降水數據降尺度研究，數據由中國氣象科學數據共享服務網站提供。此外，采用研究區內92個站點氣象站點在2000—2012年間降水觀測數據作為參照。與站點實測降水量相比，TRMM3B43數據產品普遍存在高估的現象(圖1)。已有學者研究發現，利用最小二乘方法建立TRMM3B43與站點降水數據的冪函數回歸模型，可以取得較好的校準效果[8]。

圖1 氣象站與TRMM降水數據的回歸分析圖Fig.1 Regression analysis of precipitation data from weather station and TRMM

圖1中，縱坐標表示氣象站點年總降水量，橫坐標是TRMM年均降水量，冪函數作為回歸方程時，判斷系數為0.78，擬合方程為

v=0.92·u0.95,

(1)

式中：v為氣象站年均實測降雨量；u為校準后的TRMM3B43年均降雨量。圖中虛線為y=x線，由此可見經校準模型修正后的數據較好地克服了TRMM原始數據在研究區內對降雨量高估的問題。圖2為青藏高原地區TRMM校準后的多年平均降水分布情況。

圖2 2001年青藏高原TRMM降水量空間分布圖Fig.2 Spatial distribution of TRMM precipitation in the Tibet Plateau in 2001

圖2表明，青藏高原地區降水量空間分布極其不均勻，降水量空間分布呈現自東南向西北遞減、自南向北逐漸減少的特點[14]； 局部地區受其他環境因子的影響，降水分布較為復雜。因此，考慮整個研究區內降水空間分布的整體趨勢對模型的影響，本文建立了經緯度比值指標，即

(2)

式中：ll(i,j)表示第i列第j行像元的經緯度比；lon(i,j)和lat(i,j)分別表示第i列第j行像元的經度和緯度。

2.2 歸一化差值植被指數數據

已有研究表明，NDVI與其他景觀地理學屬性，如葉面積指數[17]、地表溫度[18]、地表蒸散發[19]等有關，并且這些變量在世界不同區域均與降水有相關關系[11-12,20]。此外，甚高分辨率光譜儀(advanced very high resolution radiometer，AVHRR)的NDVI數據的時空覆蓋范圍廣，與遙感降水產品相比，空間分辨率相對較高，可以有效地提高降水數據的空間分辨率。本次研究采用美國國家航空和宇宙航行局NASA官網提供的AVHRR-NDVI日產品作為降尺度模型的輸入變量。

本文采用最大合成算法(maximum value composite，MVC)計算每個像元的月最大值，并計算每一年12個月最大NDVI數據的年平均值，應用最鄰近像元法將8 km分辨率數據重采樣為0.25°分辨率。圖3以2001年為例，青藏高原地區年平均NDVI空間分布,可知植被分布的總體趨勢，與降水量分布的總體趨勢大致相同。

圖3 2001年青藏高原NOAA-AVHRR NDVI空間分布圖Fig.3 Spatial distribution of NDVI from NOAA- AVHRR in the Tibet Plateau in 2001

2.3 SRTM DEM數據

梭雷達地形測量任務是由美國國家地理空間情報局(NGA)和NASA共同推行的國際研究項目。其雷達系統可獲得56°S～60°N近全球覆蓋范圍的數字高程模型，空間分辨率有30 m和90 m兩種。本次研究采用青藏高原地區90 m分辨率的DEM數據，并變換為0.25°尺度。數據由http: //gdem.ersdac.jspacesystem.org提供。在模型建立的過程中，除DEM外，并根據青藏高原DEM數據提取坡向和坡度數據[21]。

3 方法

3.1 降尺度方法

3.1.1 隨機森林算法

隨機森林算法(random forest,RF)在分類和預測方面具有優勢，目前已被廣泛地應用于降水和生態等諸多領域。20世紀80年代Breiman等人將分類樹方法發展成為隨機森林算法。與神經網絡算法相比，其計算量小且精度較高。該算法的優點在于： ①對多種資料，可以產生高準確度的分類器； ②可處理大量的輸入變量； ③在決定類別時，能夠評估變量的重要性。鑒于降水與多個變量存在相關關系，模型不僅需要多個輸入變量，而且有大量數據作為訓練樣本。因此，本次研究采用隨機森林算法建立降尺度模型是可行的。

相比于神經網絡和多元回歸，隨機森林算法作為集成學習方法，在分類和回歸運算方面具有優越性，該算法已經被成功用于需水預測[15]和海表面鹽度反演[16]的研究中。其建模步驟為: ①從原始訓練樣本中可重復地隨機抽取M個樣本； ②從解釋變量中不重復地隨機抽取N個樣本，從抽中的變量中選擇最能有效分割數據的變量，使分割的子集內部的變異性最小； ③依據步驟②得到的變量將數據分割為2個純度較高的子集； ④對子集重復步驟③直到分割停止，即完成了單棵分類樹的建模； ⑤重復步驟①—④X次，構建包含X棵樹的隨機森林模型； ⑥建立評價指標，檢驗模型的精度。

本文應用R軟件中隨機森林程序包作為建模工具，針對本文的研究問題進行建模，表1簡單地羅列了該程序包中包含的主要函數。

表1 隨機森林包主要函數名與功能Tab.1 Random forest package main function name and function

3.1.2 模型的假設

已有研究表明，NDVI是降水量降尺度模型的重要輸入因子。由于本文研究區內存在水體，它們對分析植被與降水的統計關系有很大的影響。因此，在建立模型的過程中不考慮NDVI小于等于零的樣本，將其視為異常值從數據中剔除。

3.1.3 模型的準備與建立

本次研究基于隨機森林算法，以2001年為例對校準后的TRMM降水數據進行降尺度方法的研究，步驟如下： ①分別將8個變量：降雨量(precipitation)，降水與歸一化差值植被指數(NDVI)、數字高程(DEM)、經緯度比值(ll)、坡向(aspect)，坡度(slope)，經度(lontitude)及緯度(lattitude)讀入數組中，使每個數組對應位置的元素值代表圖像中同一個像元的特征，并寫入同一個矩陣中； ②剔除異常值； ③將數據隨機分成大小相同的兩組，組1數據作為建模樣本，組2數據用來檢驗； ④利用R軟件的建模工具，建立Random Forest模型； ⑤輸入組2模型進行檢驗。

3.1.4 模型的應用

對上節中建立的模型進行有效性檢驗后，本節將對降水數據進行降尺度計算，具體步驟為： ①采用最鄰近像元法將NDVI，DEM，ll，slope及aspect原始數據重采樣為8 km×8 km分辨率數據； ②按照模型輸入變量的數據格式，創建經緯度和經緯度比值數據； ③剔除異常值； ④輸入已建立隨機森林，計算8 km分辨率的降水量； ⑤采用最鄰近像元法，對未進行降雨量預測的異常值像元進行降雨量插值； ⑥將上文獲得的訂正誤差與預測值求和，計算得到降尺度結果； ⑦采用最鄰近像元法獲得與氣象站位置對應像元降尺度結果，與氣象站數據進行回歸分析。

3.2 降水量的預測

降水量變化是典型的非線性時間序列，其中包含了大量的時序動態變化特征。由表2 可知，降水量與環境具有強相關性。目前常見的預測方法是將相關因子作為輸入向量，建立多元回歸預測模型，這種方法考慮了環境因子對降水時序變化的影響，但是缺乏時序動態分析，不能反映內部變化規律； 另一種常見的方法是時間序列分析，該方法充分考慮了降水變化的內部因子，卻忽略了降水與外在因子之間的關聯，進而影響預測精度。

為了提高降水量預測精度，本次研究在考慮降水變化時序動態變化特征的同時，引入了環境因子，基于隨機森林算法建立降水量的預測模型。步驟如下： 首先，需要為模型定階，即判斷降水受自身發生量影響的時效長短，采用F測驗進行逐步比較確定，本次研究采用2000—2012年均降水數據作為研究對象，預測時效期為預測時間點的前5 a； 然后，將時效時間段內的發生量作為降水內部變化的描述因子，并選擇NDVI作為外在影響因子，基于隨機森林算法建立非線性預測模型； 將2000—2004年和2005年降水數據和2005年NDVI數據作為輸入向量，建立預測模型，隨后將2001—2005年降水數據和2006年NDVI數據作為輸入向量，預測2006年降水量。

4 結果分析與檢驗

4.1 降尺度模型的檢驗

為提高訓練樣本的合理性，采用隨機方法生成樣本個數相等的訓練樣本和檢驗樣本，將表2中8個參數作為輸入量，建立隨機森林預測模型，對檢驗樣本進行預測，以2001年降水數據為例，預測結果與校準后TRMM降水數據的比較圖4所示。

圖4 隨機森林模型預測值與校準的 TRMM3B43檢驗值擬合圖Fig.4 Random forest model predictions and calibrated TRMM3B43 values

圖4中，縱坐標代表對模型輸入檢驗數據集的自變量后，輸出的降水量的預測值，橫坐標代表檢驗數據集中的TRMM校正值，兩者的判斷系數R2為 0.87。

4.2 模型的比較

鑒于賈紹鳳團隊采用線性回歸模型，對柴達木盆地的TRMM降水產品進行降尺度研究，其模型取得了較好的檢驗精度。本文建立降水量與其他自然地理變量之間的線性回歸模型，與隨機森林模型比較。如表2列出了降水與歸一化差值植被指數(NDVI)、數字高程(DEM)、經緯度比值(ll)、坡向(aspect)、坡度(slope)和經緯度(lon和lat)的線性相關性。

表2 降水與其他變量的線性相關性Tab.2 Linear correlation between precipitation and other variables

根據變量之間的線性相關系數明顯可知，降雨量與另外7個變量之間的相關性差別較大。因此，本文采用向后逐步回歸的方法，在0.25°分辨率下從模型包含所有預測變量開始，一次刪除一個變量直到會影響模型變量(Akaike information criterion, AIC)為止，即

AIC=(2K-2L)/n,

(3)

(4)

式中：k為參數的數量；L為對數似然值；n為樣本數目；sse為殘差平方和。AIC的大小取決于L和k。k取值越小，AIC越小；L取值越大，AIC值越小。k小表明模型簡潔，L大表明模型精確。因此AIC和修正的決定系數類似，在評價模型是兼顧了簡潔性和精確性。

經后向逐步回歸分析后本文建立了降水與NDVI,DEM,slope,lon,lat和ll，6個變量(X1,X2,X3,X4,X5和X6)間多元線性回歸模型(圖5)。模型預測值與檢驗值的判斷系數為0.77。鑒于模型擬合效果不如基于隨機森林的降尺度模型效果好，本文選擇后者對青藏高原地區降水數據進行降尺度研究。

圖5 多元線性模型預測值與校準的 TRMM3B43檢驗值擬合圖Fig.5 Multivariate linear model predictions and calibrated TRMM3B43 values

4.3 降尺度結果與精度檢驗

為分析降尺度結果在研究區內空間分布的準確性，本次研究將TRMM校準值(圖6)與8 km×8 km降尺度結果(圖7)進行比較分析，隨機森林輸出結果的空間分布特點呈現從南到北、從東南到西北逐漸減小的趨勢，與TRMM校準值的空間分布趨勢相同，但是在具有特殊地理環境的局部區域，例如，在珠穆朗瑪峰少雨區、祁連山脈多雨區和青藏高原腹地區域，降尺度結果不夠理想。

圖6 2001年青藏高原TRMM降水校準值空間分布圖Fig.6 Spatial distribution of TRMM precipitation calibration value of Tibet Plateau in 2001

圖7 2001年青藏高原隨機森林輸出結果空間分布圖Fig.7 Spatial distribution of random forest output in the Tibet Plateau in 2001

因此，本次研究在此建模尺度上，預測0.25°分辨率的降水量，結合TRMM校準后降水值求出模型在該尺度上的誤差分布，采用最鄰近像元方法將結果插值為8 km×8 km分辨率的誤差分布，結果如圖8所示。將8 km×8 km誤差分布與隨機森林輸出的8 km×8 km降水結果求和，即得到降尺度結果，如圖9所示。

圖8 8 km×8 km空間分辨率誤差分布圖Fig.8 8 km × 8 km spatial resolution error distribution

圖9 降尺度結果圖Fig.9 Downscaling results

本文采用回歸系數R2作為檢驗標準對降尺度結果進行驗證分析。TRMM校準值和降尺度結果與站點降水量回歸分析圖見圖10。如圖10所示，降尺度結果與氣象站點降水觀測量的判斷系數為0.89，提高了原始數據的空間分辨率。同時與0.25°校準數據(與站點數據判斷系數為0.81)相比，提高了數據與站點數據的擬合系數。因此，采用隨機森林算法對研究區進行降尺度計算，不僅將降水數據的空間分辨率從0.25°提高到8km，而且提高了降水數據的準確度。

(a) TRMM校準值與站點降水量回歸分析圖

(b) 降尺度結果與站點降水量回歸分析圖圖10 TRMM校準值和降尺度結果 與站點降水量回歸分析圖Fig.10 TRMM calibration and downscaling results with site precipitation analysis

4.4 降水量預測方法結果與檢驗

圖11分別為研究區內5個觀測站點降水量測量值和預測結果隨時間變化的曲線。5個站點依次為昌都站(31.15°N，97.17°E)、那曲站(31.48°N，92.07°E)、林芝站(29.67°N，94.33°E)、拉薩站(29.7°N，91.13°E)和日喀則站(29.25°N，88.88°E)。分析表明，在5個觀測站點處的降水量預測值有效地描述了降水的年際變化趨勢和降水量的數量級。

(a) 昌都站點 (b) 那曲站點 (c) 林芝站點

(d) 拉薩站點(e) 日喀則站點

圖115個站點觀測值與預測值的年際變化曲線

Fig.11Interannualvariationofobservedandpredictedvaluesforfivesites

采用上述預測方法，預測結果與原始TRMM降水數據的擬合系數見表3。

表3 2006—2012年預測結果與校準后的 TRMM降水量擬合系數Tab.3 Predicted results for 2006—2012 and TRMM precipitation fitting coefficient after calibration

4.5 降水數據精度對結果精度的影響

本文采用TRMM3B43產品進行降尺度研究。該產品綜合了4類相互獨立的降水數據，包括微波、近紅外等傳感器融合估算數據，以及美國國家海洋、大氣管理局和全球降水氣候中心的降水雨量計分析數據等[14]，是衛星降水數據與其他降水數據聯合反演的最佳降水產品。首先，數據產品本身存在誤差，該誤差可能由2方面引起： ①在觀測4類相互獨立降水數據的過程中引入了誤差； ②使用聯合反演算法進行數據融合的過程中引入了誤差。其次，本次研究采用站點觀測數據雖對TRMM3B43原始數據進行了校準，但所用的站點數據較少且分布不勻，導致校準模型在校準過程中引入了誤差。以上各種誤差均可影響降水數據的真實性，進而降低了降尺度結果的精度。

4.6 輸入變量對結果精度的影響

降水是受多種因素影響的氣象要素。本次研究僅僅考慮了植被和地形因素對其產生的作用，忽略了氣候帶、海陸位置、季風和人類活動等因素對降水的影響。另外，本次研究模型的輸入變量并不具有相互獨立性，在一定程度上降低了模型的有效性，因而影響了降尺度結果的準確性。

影響預測結果精度的因素主要有2個方面： ①統計模型缺乏物理機制，無法充分描述內部和外部因子對降水變化的影響。②本次研究所采用的降水數據為TRMM3B43的2000—2012年間的年均降水數據，時間尺度較短，因而影響了時間序列分析的有效性和準確性，使預測結果無法達到理想的精度。

5 結論

基于已有青藏高原的遙感降水等資料，本文根據降水與植被和地形因子的相關關系，使用TRMM3B43降水、NOAA-AVHRR NDVI和SRTM DEM等數據，采用隨機森林算法建立了0.25°尺度下的降尺度模型，并求出了該尺度下的誤差分布，采用最鄰近像元方法插值為8 km分辨率，結合8 km模型預測值計算得出了青藏高原降水分布的降尺度結果。經分析驗證，降尺度結果與地面站點降水量觀測數據的R2為0.89，高于TRMM3B43校準值與地面站點降水量觀測數據的R2(0.81)。此外，綜合考慮內部因子和外部因子對降水變化的影響，基于隨機森林算法，融合時間序列分析和多元回歸分析，建立降水量預測模型，有效地描述了降水的年際變化趨勢和降水量的數量級，并且2000—2012年均降水量預測結果與TRMM3B43降水數據擬合系數均達到0.8以上。因此本研究得出了下述認識：

1)基于下墊面因子和隨機森林算法的降尺度方法，能夠較準確地估計實際降水情況，不僅可提高衛星遙感降水產品的空間分辨率，同時可提高原始數據的反演精度。

2)融合時間序列分析的降水預測模型，可以有效地描述降水的年際變化趨勢和降水量的數量級。