利用空間隨機森林方法提升GPM衛星遙感降水質量

胡保健，李偉，陳傳法，胡占占

1.山東科技大學 測繪與空間信息學院,青島 266590;

2.湖州中核勘測規劃設計有限公司,湖州 313000

1 引 言

降水是地球能量循環與物質交換重要組成部分，是陸地生態系統重要驅動因素，在氣候變化中起著決定性作用（Markonis 等，2019；Zhou 等，2020b）。高分辨率和高精度的空間降水分布對區域水資源調度管理、災害預測防控、生態研究等領域具有重要的應用價值（Brodeur和Steinschneider，2020；Renard 等，2011；Zhang 等，2020）。然而，降水數據估算誤差在上述應用中會帶來很大的不確定性（Xie 和Xiong，2011）。如何提高降水空間分辨率和準確性是近年來生態水文等領域研究熱點和難點（Shen等，2014b）。

降水監測數據主要來自氣象觀測站、地基雷達測量、衛星遙感降水反演產品等。氣象站觀測降水可直接測得實際降水量，但受地形限制，難以反映降水的空間連續性分布（Sharifi 等，2019；Shen 等，2014a）。地基雷達觀測可以提供連續高分辨率空間降水數據，但只能反映局部區域降水，并且復雜地形區域布設困難、易受環境因素干擾多。隨著遙感技術的發展，產生了大量基于星載探測技術的衛星遙感降水產品，如TRMM、PERSIANN、GPM、CHIRPS、CMORPH（Hu等，2014；Islam等，2020；Sun 等，2016）。相對于氣象站點觀測與雷達測量降水，衛星遙感降水產品具有覆蓋范圍廣、時空連續變化效果好的特點（Zhou 等，2020a），在水文模型和生態環境研究中得到廣泛應用（Jiang等，2021；Wei等，2020）。

然而，現有衛星遙感降水產品分辨率較低（最優僅為0.05°），而且受傳感器誤差與降水反演算法不確定性等影響，現有降水產品均含有一定的系統偏差，導致其不能滿足局部尺度和生態流域尺度的應用研究（Immerzeel 等，2009）。因此，國內外研究者們相繼提出了多種點面融合方法，以借助地面觀測站降水對衛星降水產品修正，如最優插值法（OI）（盧新玉 等，2017；潘旸 等，2012）、貝葉斯模型（Fu等，2016）、地理差異分析（GDA）（Jongjin等，2016）、地理加權回歸（GWR）（Lu等，2019）、K最近鄰模型（KNN）（王玉丹 等，2016）等。這些方法雖然有效提高了衛星降水產品的精度，但仍有一些問題需要解決：（1）融合后的降水產品分辨率與原有的衛星降水產品相同，不利于局部地區水文模型研究；（2）一些方法只考慮了實測降水背景場，未考慮與降水有關的環境因素；（3）在進行偏差校正時，粗分辨率的衛星降水產品與基于點的觀測站實測降水存在尺度差異，導致點面融合時捕捉不到局部降水細節信息。因此，在點面融合之前對遙感降水產品進行空間降尺度是十分必要的。

近年來，國內外學者們提出和發展了多種遙感降水產品降尺度方法，如指數回歸（ER）（Immerzeel等，2009）、多元線性回歸（MLR）（嵇濤 等，2015；Jia 等，2011；馬金輝 等，2013）、GWR 方法等（胡實 等，2020）。隨著研究深入，一些研究者們開始嘗試使用機器學習方法對遙感降水進行降尺度研究。大量研究表明，隨機森林RF（Random Forest）模型在處理多維數據中具有良好的表現，且能夠很好的捕捉自變量與因變量之間非線性關系（Baez-Villanueva 等，2020；Njuki 等，2020），具有良好的應用前景（Jing 等，2016；Ma 等，2018）。然而，上述降尺度研究均將RF 模型作為一種簡單的統計工具，忽略了相鄰遙感（站點）降水數據的空間相關性。

基于以上討論，本文提出了一種顧及空間相關性的雙階段衛星降水產品質量提升方法，即降尺度與點面融合相結合的思路，以提升衛星降水產品的空間分辨率與精度。該方法以RF 為基礎模型，構造了顧及空間相關性的空間隨機森林模型SRF（Spatial Random Forest）。與現有方法相比，新方法優勢包括：（1）考慮了降水之間空間相關性；（2）降尺度與點面融合兩階段都采用了SRF模型，且均融入了環境變量影響。以四川省為研究區域，選取最新GPM IMERG V06B 遙感降水數據作為數據源，使用SRF-SRF 方法對GPM 降水數據進行降尺度與點面融合處理，并將計算結果與經典方法比較，驗證新方法的可行性和高效性。

2 研究區域和數據源

2.1 研究區概況

四川省位于中國西南部（97°21′E—108°31′E，26°03′N—34°19′N），總面積約為4.86×105km2，介于青藏高原與長江中下游，是兩者的過渡地帶。研究區域地形地貌復雜多樣，涵蓋了山脈、高原、丘陵、平原和盆地等地形，地形起伏走勢西高東低。受地形因素和空間位置影響，該區域降水夏季多冬季少，其中夏季降水占全年80%—90%，且受川西垂直地勢起伏影響，四川盆地降水充沛，年降水量800—1600 mm，呈現由東到西降水逐漸減少的變化特點（Yang 等，2020；Lai 和Gong，2017）。圖1為研究區域地形及氣象站點分布。

圖1 四川省氣象站分布Fig.1 Distribution of weather stations in Sichuan Province

2.2 數據源

（1）氣象站數據。氣象站觀測數據由中國氣象局氣象數據中心（http://data.cma.cn/［2021-01-01］）提供。選取四川省內2015 年—2019 年期間156 個逐日氣象站點降水數據（圖1），所有數據經過內部一致性檢查、極值檢查以及空間一致性檢查等嚴格質量控制后（Song 等，2004），經累積計算得到每月降水量。

（2）衛星降水數據。GPM 降水數據從美國宇航局地球科學數據中心（https://pmm.nasa.gov/［2021-01-01］）獲取，是美國NASA（美國國家航空航天局）和日本JAXA（日本宇宙航空研發機構）共同研發的繼TRMM 之后全球降水測量衛星，搭載全球先進的DPR（Dual-frequency Precipitation Radar）和GPM GMI（GPM Microwave Imager）系統，對微小降雨探測能力比前一代產品TRMM 更精準，覆蓋范圍（60°S—60°N）更廣闊。本文選用月尺度0.1°×0.1°分辨率的降水產品GPM IMERG V06B（簡稱GPM）數據集作為本研究所采用的數據源，時間跨度為2015年1月—2019年12月。

（3）DEM、NDVI 以及地表溫度（LST）數據。研究采用的DEM 數據為SRTM（Shuttle Radar Topography Mission）DEM V4.1 數據集，從地理空間數據云（http://www.gscloud.cn/［2021-01-01］）獲取，空間分辨率為90 m，經像素平均重采樣生成1 km 分辨率，然后在Arcgis 中利用表面分析工具獲得坡度、坡向、地形起伏度數據。NDVI（MOD13A 3）、LST（MOD11A2）數據從美國航空局（https://ladsweb.modaps.eosdis.nasa.gov/［2021-01-01］）中獲取，空間分辨率為1 km，NDVI 時間分辨率為月，LST時間分辨率為8 d。其中，LST數據包括白天地表溫度（LSTd）與夜間地表溫度（LSTn）以及白天與夜間的溫度差（LSTd-n），由每8 天LST 求平均得出月平均溫度LST。

3 研究方法

本文構建了一種雙階段衛星降水產品質量提升方法。第一階段，通過考慮降水與其他因素（如地形、NDVI、地表溫度、經緯度等）之間的關系，借助SRF 將10 km GPM 衛星降水產品降尺度到1 km 空間分辨率（D_GPM）。第二階段，將實測站點降水數據與D_GPM 及上述自變量（如地形、NDVI、地表溫度、經緯度等）相結合，再次借助SRF 生成高精度高空間分辨率的降水數據。該方法的流程圖如圖2所示。

圖2 雙階段降尺度流程圖Fig.2 Two-stage downscaling flowchart

3.1 隨機森林（RF）模型

隨機森林（RF）算法是Breiman 于2001 年提出的一種基于bagging 的集成學習方法，通過構造多種決策樹來處理自變量與因變量之間的關系，可用作數據的分類和回歸預測。RF 通過構建大量的樹模型，對多種特征值的重要性進行整合篩選，充分考慮不同特征值之間的重要性，選擇最優樣本特征值進而找到最優解，并求得所有預測值的平均值作為最終預估值。與傳統回歸預測方法相比，RF 可以處理復雜多維的特征值，回歸預測更準確，采取樹狀重復放回抽樣，避免過擬合，且不必考慮線性回歸中的多重共線性問題，穩健性更好，在多種領域得到廣泛應用。RF 模型的通用公式如下表示：

式中，P（s0）為s0處的降水預測值，Xi（s0）（i=1，2，3，…，k）為s0處的自變量，k為自變量個數，ε為預測誤差。

3.2 顧及空間相關性的隨機森林（SRF）模型

為充分考慮鄰近降水之間的空間相關性，本文構建了顧及空間相關性的空間隨機森林（SRF）模型。其中，SRF 是在RF 的基礎上引入克里金估計值作為輸入變量，其通用公式如下表示：

式中，Ps（s0）為s0處的克里金降水估計值。

克里金插值法是一種無偏最優估計值方法（史文嬌 等，2012；Chen 和Li，2019；Kim 等，2013），對于預估點s0處的降水值PS（s0）可通過搜索該點周圍鄰近n個站點的降水P（x）經線性加權求得，其計算公式為

式中，PS（s0）為s0處的降水估計值，λi為克里金權重系數，n為周圍鄰近站點數，P（si）為站點si處的降水值。

克里金權重系數λi不僅取決于周圍鄰近降水站點與預估點的距離，還取決于鄰近點空間分布，其值是通過求解克里金目標函數實現，即PS（s0）對真實值P（s0）的估計為無偏（式（4））且方差最?。ㄊ剑?））。

由于區域性降水分布具有較大的空間相關性（Sekuli? 等，2020），因此，借助克里金插值可充分考慮樣本之間的空間相關性。

3.3 SRF雙階段降尺度

本研究的雙階段衛星降水產品質量提升方法大致分為SRF 降尺度與SRF 點面融合兩部分，具體步驟如下：

（1）克里金插值。將10 km 分辨率的GPM 數據通過克里金插值獲得10 km 分辨率的GPM10kmkri和1 km 分辨率的GPM1kmkri，對站點實測降水數據（RGS）進行克里金插值獲得1 km分辨率RGS1kmkri。

（2）將1 km 空間分辨率的所有自變量X1km（包括NDVI、LSTd、LSTn、LSTd-n、DEM、坡度、坡向、地形起伏度以及經緯度）分別使用像元平均重采樣至10 km空間分辨率，記為X10km。

（3）以步驟（2）中重采樣自變量X10km以及步驟（1）中GPM10kmkri為自變量，原始GPM 為因變量建立SRF降尺度模型：

式中，PGPM10km(s0)為s0處的GPM降水值。

（4）將所有1 km 空間分辨率的自變量X1km輸入到SRF 降尺度模型獲得1 km 空間分辨率的衛星降水D_GPMsrf：

（5）對降尺度數據進行點面融合，即建立以RGS為因變量的校正模型，自變量包括以步驟（4）中的降尺度結果D_GPMsrf、1 km空間分辨率下的自變量X1km以及RGS1kmkri，即：

式中，RGS(s0)為站點s0處的實測降水。

（6）SRF 模型訓練以后，將全部1 km 分辨率自變量輸入到SRF 融合模型中，獲取融合校正后的1 km降水數據DC_GPMsrf。

3.4 評價方法

為評價新方法（SRF-SRF）的可行性和高效性，本文選取了7 種方法與其進行比較，包括3 種傳統方法，即GWR、RF、BPNN；3 種基于SRF 框架下的方法，即GPM 經雙線性插值降尺度后使用SRF方法進行點面融合方法（Bi-SRF）、對GPM 使用SRF 降尺度后利用站點實測降水進行地理差異分析校準方法（SRF-GDA）、對GPM 年降水使用SRF 降尺度后按月比例分解并經SRF 校正方法（SRFdis）；最后一種方法為借助克里金直接對站點降水插值（記為Kriging）。

本文采用十折交叉驗證方法驗證所有方法的性能。選取的精度評價指標包括平均絕對誤差（MAE）、均方根誤差（RMSE）、相關系數（CC）（Lu等，2020）。3種評價指標公式如下：

式中，Poi、Psi分別代表第i個站點降水實測值和該站點對應的降水預估數據，分別代表站點實測降水的算術平均值和該站點對應降水預估數據算術平均值，n表示氣象站點個數。相關系數CC表示兩種數據相關性，值越接近1 表示相關程度越好。RMSE 和MAE 用于評價預估降水與實測數據之間的誤差，值越小表示精度越高。

4 結果與分析

采用SRF-SRF 方法對四川省GPM 月降水數據進行降尺度與點面融合研究，并將該方法的結果與經典方法的結果在月、季、年3種時間尺度上進行比較。

4.1 月尺度分析

表1 顯示了各模型在2015 年—2019 年全部月份上的計算精度。結果顯示，在所有模型結果中，GWR 與BPNN 效果最差，可能是這兩種方法沒有有效處理降水與環境變量之間的復雜關系?；谡军c的kriging比RF方法好，且所有考慮空間相關性模型（包括SRF-SRF、Bi-SRF、SRF-GDA、SRFdis）計算結果精度都優于傳統的方法，說明了引入降水之間的空間相關性信息對模型精度提升尤為重要。與原始GPM相比，SRF-SRF的MAE和RMSE分別降低了19.51%、16.35%。整體而言，SRF-SRF 精度優于Bi-SRF、SRF-GDA、SRFdis方法。

表1 2015年—2019年全部月份各模型精度對比Table 1 Precision comparison of each model in all months from 2015 to 2019

各個方法在月尺度上的MAE、RMSE、CC 箱線圖如圖3 所示。結果顯示，BPNN 方法在所有模型中效果最差，MAE、RMSE、CC 等3種指標的中位數分別為22.66 mm、30.48 mm、0.64，其次為GWR、RF、Kriging 方法，這與全部月份精度對比（表1）一致?；赟RF的4種方法結果精度都優于傳統方法，其中SRF-SRF方法的MAE、RMSE、CC等3 種指標的中位數分別為15.66 mm、21.03 mm、0.81，而其他3 種基于SRF 方法的3 個指標分別為15.83—16.15 mm、21.41—22.27 mm、0.77～0.79，進一步說明了SRF-SRF的高效性。

圖3 2015年—2019年月尺度精度指標箱線圖Fig.3 Box plot of monthly scale accuracy indicators from 2015 to 2019

圖4 為GWR、BPNN、RF、Kriging、SRFdis、SRF-SRF 在月尺度上各站點RMSE 分布圖。由于受地勢和氣候影響（圖1），來自印度洋水汽在四川中部盆地與高原相交地段受阻凝結形成降水，故雅安與樂山地區降水偏多，導致其較大的RMSE。GWR、BPNN、RF 方法在各個站點的RMSE 均高于基于SRF 的模型，尤其在四川盆地降水較多地帶。這是因為在盆地區域站點分布較為均勻，考慮站點之間的相關性可以有效校正該區域的誤差。由圖4（e）和 圖4（f）表明，與SRFdis 相比，SRF-SRF 方法在四川東北地區與四川盆中地區略好一些。

圖4 2015年—2019年月尺度RMSE站點分布圖Fig.4 RMSE distribution of all sites on a monthly scale from 2015 to 2019

4.2 季尺度分析

表2 顯示了各模型在不同季節上的精度表現。四川地區氣候普遍表現為冬暖夏熱，夏季降水多，春秋冬降水少，所以在夏季降水的GPM 反演降水誤差較高。具體而言，GWR 與BPNN 方法在冬季的表現要優于其他季節，可能該地區冬季氣候較為溫暖，降水與不同環境因子之間相關性較強。3種傳統方法中，RF優于GWR、BPNN，而所有基于SRF 方法（除冬季的SRF-GDA 方法外）的精度都優于這3 種方法。SRF-GDA 在冬季的效果較差，主要是因為冬季降水較少，且在進行點面融合時，未考慮站點以外的環境變量。SRF-SRF 方法在這所有模型中方法效果最好，其中在冬季提升效果明顯；相較于原始GPM，SRF-SRF 的MAE 降低了40.69%，RMSE降低了44.24%，CC提高了43.36%。

表2 2015年—2019年季尺度各模型精度對比Table 2 Accuracy comparison of various models on a seasonal scale from 2015 to 2019

4.3 年尺度分析

圖5顯示了各模型在年尺度上的精度表現。由于四川省從2015 年—2018 年年降水量逐年增加，且在2018 年降水量最多，所以各方法在2018 年的MAE、RMSE 最大。所有方法中，BPNN 最差，其次為GWR，其與月尺度和季尺度表現效果一致。與Kriging 與GWR、BPNN 相比，RF 略好一些；基于SRF 模型效果表現均優于傳統的方法。整體而言，SRF-SRF 精度最優，而SRF-GDA 方法在所有基于SRF 模型中表現最差，這主要是因為后者在降水校正時未考慮環境變量有關。

圖5 2015年—2019年年尺度精度指標Fig.5 2015—2019 annual scale accuracy indicators

4.4 降水的空間分布特征

圖6 顯示了各方法在2018 年7 月的降水空間分布圖，其中，該月份為5 年內降水最多的月份。從降水分布可以看出，四川省地區降水分布空間差異明顯，與地形走勢較為接近，即降水多分布于四川中東部地勢較低地區，而川西地區由于地勢突起降水相對較少。原始GPM（圖6（a））降水有效捕捉了川西地區的降水空間分布與降水量，而在四川中部地區卻低估了降水，這可能是因為四川中部地形復雜，遙感衛星監測時受四川中部地勢和潮濕氣候影響，導致反演的衛星降水與實測降水具有一定的誤差。經降尺度融合后，各種方法明顯改善了中部地區的低估情況，空間分布整體上呈一致性。就空間分辨率而言，原始GPM（圖6（a））降水數據分辨率較為粗糙，克里金插值站點結果（圖6（b））雖然分辨率提升，但圖像比較模糊，細節特征不明顯；其他方法降水空間分布細節信息豐富，空間分辨率得到顯著提升。GWR 方法（圖6（c））的降水空間分布存在較多的噪聲，即在部分區域呈現斷崖式降低，不符合降水連續性分布的特點；RF 方法（圖6（d））的降水在四川東部與四川南部呈現部分塊狀分布，但整體上相對于GWR 更平滑一些；SRFSRF 方法（圖6（e））更能刻畫出四川中部地區降水的空間分布情況，且減少了RF 預測模型結果的塊狀分布，降水分布在空間上更具有連續性，符合降水分布的特點。

圖6 2018年7月降水空間分布圖Fig.6 The spatial distribution of precipitation in July 2018

5 討 論

5.1 RF模型變量重要性分析

降水是大氣與環境相互作用的結果，因此降水與當地的地形、植被等因素存在著較高的相關性。地理空間位置與地形因素以及植被因素一直以來是研究降水降尺度的常用變量（李凈和張曉，2015），然而，在一些高海拔地區，受積雪覆蓋等影響，NDVI 和地形因素與降水之間的關系往往不太明顯。LST在植被稀少的地區與降水之間存在著顯著關系（Jing 等，2016），即降水能影響局部地區的溫度變化。因此本方法考慮了空間位置、地形因素、植被因素與地表溫度作為降水降尺度的環境變量?；赗F 模型的相對變量重要性（圖7）分析表明，Kriging 插值降水的重要性最大，其次為SRF 降尺度值，說明研究區域站點實測降水之間具有顯著的相關性。地表溫度的重要性要高于NDVI，可能研究區域中高海拔地區較多，溫度較低不利于植被生長，難以有效反映降水和NDVI 關系。地形因子中坡向的重要性最低，其次是坡度。而空間位置對降水的影響也很重要，這與研究區域所處的氣候有很大關系。

圖7 SRF模型中各變量的重要性Fig.7 The importance of variables in the SRF model

5.2 NDVI滯后性影響分析

有研究表明，NDVI 可能對降水存在3 個月以上的滯后性（胡實 等，2020；Karbalaye Ghorbanpour等，2021）。為驗證這種長時間滯后性是否對降尺度融合產生影響，本研究將年均NDVI 融入到年降水降尺度中，然后按月降水占的比例獲得月降水降尺度結果，最終進行SRF 點面融合（即SRFdis）。圖8給出了2015 年—2019年月尺 度SRF-SRF 降水預估值與SRFdis 降水預估值的散點分布圖。兩種方法的離散程度與實測降水的擬合程度都接近于1∶1 線，與原始GPM（表1）相比均顯著提高了GPM 降水的準確性。SRF-SRF 與SRFdis 相比，前者的MAE 和RMSE 值比后者分別降低了0.36 mm 和0.76 mm。因此，與年平均NDVI作為環境變量相比，選擇月NDVI 值可以有效提高計算結果精度。

圖8 SRF-SRF與SRFdis降水預估值與實測降水散點分布圖Fig.8 Comparison between SRF-SRF and SRFdis

6 結 論

為提高GPM IMERG 遙感降水產品空間分辨率與準確性，本文以RF 模型為基礎，構造了一種顧及空間相關性的空間隨機森林模型，主要包括降尺度和點面融合兩個階段。首先使用SRF 模型并結合降水之間的空間相關性以及環境變量對GPM IMERG 降尺度，然后以實測降水產品為基準再次基于SRF 模型對降尺度后的結果點面融合，獲得高空間分辨率和高精度降水產品。實驗結果表明：

（1）在復雜地形區域考慮降水空間相關性可有效提升降水產品空間分布細節特征和降水的準確性，并且通過對RF 模型變量重要性分析得出kriging 插值信息重要性最大，說明了研究區域內降水之間的空間相關性不可忽視。該方法解決了研究區域內因降水產品空間分辨率粗糙和準確度低等導致難以開展精細化水文研究等難題，為降水產品精細化研究提供了技術支持。

（2）為驗證本文方法的有效性，將該方法實驗結果與其他經典方法結果加以對比，結果表明不同時間尺度上該方法表現效果均優于其他方法，并且本文方法的預估值與實測降水更具一致性，說明在降尺度與點面融合兩階段都考慮與降水相關的影響因素可提高結果精度。

后續研究中，將使用其他高時空分辨率降水驅動因素（如土壤濕度、地表溫度、風速等）并融合不同遙感降水產品，進一步提升遙感降水產品時空分辨率（如逐日、逐小時）和精度。