基于GWR模型的長江流域TRMM數據降尺度

張寒博, 徐 勇, 竇世卿, 靖娟利, 張 楠, 張偉東

(1.桂林理工大學 測繪地理信息學院, 廣西 桂林 541006; 2.黑龍江省農墾科學院, 哈爾濱 150038)

降水是地球水循環中重要的一環，在不同時空尺度的大氣過程中也扮演著不可或缺的角色，是氣象學、生態學和水文學的重要參數，推動全球物質循環和能量交換。傳統的區域降水數據一般是基于氣象站觀測值通過以點帶面的空間插值方法獲取[1-3]，因此需要站點數據滿足分布均勻，且分布密度較高[4]，實際較難達到要求。

近年來，衛星遙感技術發展迅速，通過遙感數據進行降水探測已成為空間化降水數據的重要來源[5]。1997年TRMM(Tropical Rainfall Measuring Mission)衛星的成功發射，更是為與降水有關的研究提供了更豐富的數據源[6]。TRMM能提供較為精確的高時間分辨率降水數據，反映降水的時間變化特征及空間分布特征[7-8]。國內外學者在不同時空尺度、不同地區對TRMM數據精度已進行了檢驗，研究表明TRMM數據與氣象站點數據、雷達觀測數據有較好的一致性[9]。谷黃河等[10]通過TRMM數據與氣象站實測數據在時間和空間分布上的對比分析，得出了TRMM數據在長江流域具有較高精度。

在各類遙感降水產品中，TRMM數據的空間分辨率已較為出色，但仍未滿足區域空間精度要求，因此，對TRMM數據進行降尺度研究十分必要。通過降尺度技術實現遙感數據空間分辨率的提高，有助于推進中小尺度區域降水時空變化特征的研究。目前，國內外學者通過不同方法獲得降尺度TRMM數據。常用的降尺度模型有兩種：一般線性回歸模型和地理加權回歸模型，一般線性回歸模型采用經典統計學理論中的最小二乘法(Ordinary Least Squares，OLS)，具有完善的統計推斷方法，但由于降水與地表特征之間存在空間非平穩性和尺度依賴，因此忽略了數據的局部特性[11-13]。地理加權回歸模型不僅強調了空間局部特征，還能反映降水的空間異質性[14]。Chen[15]基于GWR模型對華北地區TRMM 3B43數據進行了降尺度，結果表明通過GWR模型可以得到高分辨率降水數據；姬世保等[16]選取NDVI、坡度、坡向等作為自變量，降水作為因變量建立GWR模型和多元線性回歸模型，將黔桂喀斯特山地TRMM數據從0.25°×0.25°降尺度為0.05°×0.05°，并通過交叉驗證法對降尺度結果進行優選。研究結果表明GWR模型精度優于多元線性回歸模型；李瓊等[17]以黃河源區降水數據為因變量，經緯度、地形因子主成分變量為自變量，采用逐步回歸、BP神經網絡、GWR地理加權回歸3種方法對TRMM 3B43數據進行降尺度，研究表明GWR地理加權回歸降尺度效果最好；Fung等[18]通過建立EVI值和降水量之間的關系進行線性回歸，研究表明EVI和降水量的關系比NDVI更為密切；目前國內對遙感降水數據的降尺度研究主要是通過建立單要素或多要素與遙感降水數據的全局回歸模型為主[19]，且大多以NDVI為解釋變量。

本文在考慮地形因子的同時，通過不同植被指數反演TRMM降尺度數據，分別以EVI，NDVI為解釋變量，DEM、坡向為控制變量，基于GWR模型對長江流域TRMM 3B43數據進行降尺度研究。并將氣象站觀測降水量作為驗證數據集，分別對原始TRMM數據、降尺度數據(表1)進行精度驗證與結果與分析。

表1 降尺度數據信息

1 研究區與數據

1.1 研究區概況

長江流域(圖1)位于24°27′—54′N、90°33′—122°19′E，流域面積180萬km2，約占我國國土面積的19%，發源于青藏高原格拉丹東雪山西南側，東臨太平洋。干流流經11個省份，流域形狀表現為東西長、南北窄，地勢西高東低，起伏較大，呈三級階梯狀[ 10]。長江流域水資源豐富，是中國水資源配置的戰略水源地，降水主要受東南季風和西南季風影響，在時間上具有明顯的年際和季節特征，在空間上呈現出顯著的空間異質性[20]。此外，長江流域不僅在區域水循環、氣候變化等方面研究上有著重要作用，同時還是中國經濟發展的主驅動軸[21]。

圖1 研究區及氣象站點分布

1.2 數據來源與預處理

本文采用2001—2019年的TRMM 3B43(Version 7) 作為降水數據源，其空間分辨率為0.25°×0.25°(約27.5 km×27.5 km)，時間分辨率為1個月，空間范圍為(50°S—50°N，0°—180°—0°)，格式為HDF。NDVI，EVI數據來自terra衛星MODIS傳感器的植被指數數據MOD13A3，空間分辨率為1 km×1 km，時間分辨率為1個月。3種數據都通過NASA數據庫免費獲得(https:∥search.earthdata.nasa.gov/)。DEM數據從地理空間數據云(http:∥www.gscloud.cn/)SRTM數據集獲取，空間分辨率為90 m×90 m。降水量地面觀測數據選取長江流域2001—2019年147個地面站(圖1)氣象數據日值數據集，數據內容包括各站點經緯度、海拔、日降水量等。數據來源于中國氣象局國家氣象中心(http:∥cdc.cma.gov.cn/home.do)。

通過MRT(MODIS Reprojection Tool)在Java環境下對NDVI，EVI數據進行批量波段提取、拼接、投影變換、區域裁剪、無效值剔除、最大值合成等操作，得到2001—2019年NDVI月數據、EVI月數據。并通過ArcGIS對一年12期的NDVI，EVI數據求平均，分別得到年均NDVI、年均EVI數據；TRMM數據借助ENVI軟件進行批量投影轉換、旋轉、格式轉換、單位轉換等預處理，將一年12期的TRMM數據進行求和，獲得TRMM年累計數據；DEM數據進行投影轉換，再采用最鄰近內插法將其重采樣為1 km×1 km，并利用ArcGIS計算出坡向數據；降水量地面觀測數據經質量控制，剔除無效數據，對缺失值進行插補后累加得到各站點月、年時間尺度的降水總量。

2 研究方法

2.1 GWR模型介紹(原理)

2.1.1 地理加權回歸模型GWR 地理加權回歸模型(Geographically Weighted Regression，GWR)最早是由Brunsdon等[22]于1996年提出的一種用于量化空間異質性的局部參數估計方法[23]，是傳統回歸模型的擴展。地理加權回歸模型應用了局部回歸的思想，遵循地理學第一定律，通過引入數據的空間位置計算數據在局部回歸方程中的權重。在考慮相鄰點的空間權重下，通過估算每一位置的因變量與解釋變量的參數來建立回歸模型[24-25]。其基本公式如下：

(i=1，2，…，m)

(1)

式中:yi為第i個樣本點的降水量,作為因變量,共有m個點;xit為第t個解釋變量的第i個樣本點的觀測值;(ui,vi)表示第i個樣本點的經緯度坐標;β0(ui,vi)為第i個樣本點的常數項回歸參數;βt(ui,vi)為第t個影響因子對第i個樣本點的線性回歸參數；ε(ui,vi)為模型在第i個樣本點所計算出來的殘差值。

2.1.2 高斯(Gauss)函數法原理 空間權重矩陣的計算是GWR模型求解的關鍵。高斯函數法可以較高的精度計算出權重比，并基于地理學第一定律，通過距離和權重呈現反比關系，構建出一個單調遞減的連續性函數式，公式如下：

(2)

式中：Dit為研究區點與點之間的距離;b為帶寬,用來描述權重與距離函數關系式的非負衰減參數; exp[-(Dxi/b)2]為以e為底的指數函數。

2.1.3 CV法則原理 帶寬GWR的大小對GWR模型結果有較大影響，國際上最普遍的方法就是用Cleveland和Bowman提出的交叉驗證(cross-validation，CV)方法來確定，公式如下：

(3)

2.2 降尺度方法

本文利用TRMM數據和NDVI，EVI，DEM，坡向影響因子基于MATLAB環境建立函數關系，來實現TRMM數據的降尺度，主要步驟如下：

(1) 將1 km×1 km分辨率的NDVI，EVI，DEM，坡向數據重采樣為0.25°×0.25°。將數據分為高分辨率組數據(1 km)和低分辨率組數據(0.25°)。

(2) 將統一尺度的低分辨率組數據中的TRMM數據(0.25°)為因變量，NDVI數據(0.25°)/EVI數據(0.25°)為解釋變量、DEM數據(0.25°)、坡向數據(0.25°)為控制變量，建立GWR訓練模型。從回歸模型中得到常數項，解釋變量、各控制變量對應系數及殘差結果。

(3) 將上步的模型結果柵格化，并將常數項，解釋變量、控制變量系數重采樣為1 km，對殘差結果通過反距離權重法插值得到高分辨率的殘差(1 km)。

(4) 將高分辨率組數據按照GWR模型原理進行回代，將解釋變量系數、各控制變量系數(1 km)與同分辨率的解釋變量、控制變量相乘，并與常數項(1 km)相加，得到預測降水數據(1 km)。

(5) 將預測降水數據與同分辨率殘差數據相加得到最終的降尺度降水數據(1 km)。

2.3 插值方法的選取

通過GWR模型得到的預測數據與原始TRMM數據存在著一定的殘差，需要借助插值方法對殘差數據進行尺度轉換處理。本文選取研究區域5 367個樣點，分為3 634個訓練點和1 733個驗證點，通過交叉驗證對自然鄰點法、樣條函數法、普通克里金法、反距離權重法4種差值方法進行比較，并以均方根誤差(RMSE)、決定系數(R2)作為驗證指標。交叉驗證結果見表2，本研究選擇反距離權重插值法進行殘差插值。

表2 交叉驗證結果

2.4 評價指標

本文以氣象站實測數據為“真實值”，引入決定系數(R2)、均方根誤差(Root Mean Square Error,RMSE)、相對誤差(BIAS)及平均絕對偏差(MAE)4個指標對降尺度結果進行精度驗證。R2評定實測值與降尺度降水數據的相關性，其值在0到1之間，值越大相關性越高；BIAS反映“真實值”與降尺度降水數據的偏離程度，值越接近0精度越高；RMSE評定誤差的整體水平；MAE則評定誤差的實際情況。公式如下：

(4)

(5)

(6)

(7)

3 結果與分析

3.1 TRMM數據適用性分析

衛星降水數據與氣象站點實測數據在不同時空尺度存在偏差。為了保證TRMM 3B43數據降尺度的可行性與降尺度結果的可靠性，首先對TRMM數據的適用性進行分析。以長江流域2001—2019年147個氣象站數據為自變量，對應年尺度TRMM降水數據為因變量，進行一元線性回歸分析。并基于氣象站點通過建立Thiessen多邊形對研究區進行分割，得到站點相關系數的空間分布。由圖2可知，長江流域TRMM降水數據與實測降水數據之間的相關系數為0.46～0.97，且大多數站點在0.75以上，東部相關系數整體優于西部。通過相關系數顯著性檢驗，所有站點通過了95%顯著性檢驗，式中135個氣象站在p<0.01水平下顯著。這說明兩種數據之間具有較好的線性相關性，在長江流域進行年尺度降水的降尺度研究具有一定的科學合理性。

圖2 氣象站點相關系數分布

3.2 降尺度結果與驗證

3.2.1 多年平均降水降尺度結果 對2001—2019年的TRMM 3B43年累計降水數據、年均NDVI數據、年均EVI數據分別求平均得到多年平均數據，通過GWR模型獲得降尺度TRMM數據。由圖3可知：TRMM原始數據年均降水范圍為251～2 145 mm，NDVI_降尺度TRMM數據范圍為60～2 415 mm，EVI_降尺度TRMM數據為81～2 599 mm；洞庭湖支流區域降尺度前后細節表征更突出。整體空間分布特征上，TRMM數據在降尺度前后趨于一致，呈現由東向西遞減的特點，其中東南地區年降水量一般大于1 600 mm，局部地區大于2 000 mm。原始TRMM數據空間分辨率相對較粗糙，局部特征不明顯，經過GWR降尺度計算，空間分辨率從0.25°提升至1 km，其降水范圍大于原始TRMM數據，細節性增強，能更好體現降水特征。

圖3 2001-2019年平均年降水量空間分布

為進一步驗證TRMM數據降尺度結果的準確度和精度，將研究區內147個氣象站點多年平均觀測數據分別與對應的降尺度前后TRMM數據進行一元線性回歸分析(圖4)，以R2，BIAS，RMSE和MAE共4種統計指標對降尺度結果進行驗證。驗證結果表明：TRMM降水數據與氣象站點觀測數據之間的R2，BIAS，RMSE以及MAE分別為0.911 6，0.061 56，131.862，102.557，TRMM降水數據與氣象站點觀測數據之間具有明顯的線性相關性，且原始TRMM數據的降水量整體比站點觀測數據降水量偏大。

與TRMM原始數據的降水量相比，兩種降尺度數據的R2，BIAS均有所下降，RMSE和MAE略微上升。NDVI_降尺度降水數據R2指標值為0.906 2，更接近于0.911 6；BIAS值增加了0.019 6，RMSE指標值提升了3.371 mm，MAE提升0.51 mm；EVI_降尺度降水數據R2指標值為0.900 7，減少了0.010 9，BIAS值增加0.036 5，RMSE指標值提升了7.059 mm，MAE提升4.029 mm。整體來看，NDVI_降尺度模型相比EVI_降尺度模型精度更高。

圖4 降尺度前后TRMM數據與實測降水數據散點圖

3.2.2 年尺度降水降尺度 把2001—2019年氣象站點的實測降水數據和對應的原始TRMM 3B43降水數據、NDVI_降尺度TRMM數據、EVI_降尺度TRMM數據進行逐年檢驗和對比分析。圖5為3種數據與氣象站點實測降水數據精度檢驗得出的4種參數對比圖。由圖5A可知：2001—2019年3種降水數據各年R2曲線趨勢表現為明顯的一致性，19 a間原始TRMM數據與氣象站實測降水量的R2普遍相對較好；NDVI_降尺度降水量與氣象站實測降水量的R2與前者幾乎接近，2010年、2019年比前者相關性更好；EVI_降尺度降水量與氣象站實測降水量的R2相比前兩種表現稍差。降尺度前后數據的R2有所下降，其原因是在尺度轉換過程中會產生誤差，導致降水量估計有偏差。圖5B-5D為3種數據與站點實測數據的BIAS，RMSE及MAE值的對比，其中BIAS原始TRMM數據表現最好，NDVI_降尺度降水數據次之，EVI_降尺度降水數據相對較差；2006年、2010年、2012年、2015年、2019年NDVI_降尺度降水數據與氣象站點實測數據的RMSE和MAE最小，能更好反映真實降水，其余各年略微低于原始TRMM數據；各年內EVI_降尺度前后TRMM降水數據與氣象站實測數據間的RMSE和MAE都大于另兩種數據，精度相對較差。整體而言，DVI_降尺度結果優于EVI_降尺度結果，在提高空間分辨率的同時較好地保持數據的精度，能更真實地反映研究區域內的降水特征。

圖5 2001-2019年精度檢驗結果

3.2.3 典型年降水降尺度結果與驗證 2001—2019年年均降水量為1 109.31 mm，2011年降水量924.91 mm為近19 a最低，2016年降水1 287.25 mm為近19 a最高。以2011年作為干旱年、2016年作為濕潤年，進行典型年份降水降尺度結果與分析。按上述降尺度步驟對2011年、2016年TRMM數據進行降尺度，得到相應年份的NDVI_降尺度預測TRMM和EVI_降尺度預測TRMM(圖6B)，并對殘差插值得到(圖6C)1 km殘差值，最后得出相應年份的NDVI_降尺度TRMM數據和EVI_降尺度TRMM數據(圖6D)。由圖6可知，兩種降尺度結果空間分辨率都有大幅度提升，細節性增強，能更好地體現降水特征。典型年份(2011年、2016年)的降水量空間分布趨勢與多年年均降水相近，整體呈現從東向西逐漸減小的趨勢，其中，干旱年份降水主要分布在長江流域中北部和東北部，濕潤年份降水主要分布在長江流域東南部。2011年原始TRMM年均降水范圍為235～1 874 mm，NDVI_降尺度結果降水范圍為153～1 964 m，EVI_降尺度結果降水范圍為145～2 001 mm；2016年原始TRMM年均降水范圍為183～3 022 mm，NDVI_降尺度結果降水范圍為131～3 308 mm，EVI_降尺度結果降水范圍為116～3 315 mm。兩種降尺度結果都使降水區間范圍擴大，NDVI_降尺度數據范圍變化相對較小，更加符合實際降水特征。

利用長江流域2011年、2016年147個氣象站點實測數據對兩種降尺度結果進行驗證。從表3可以看出濕潤年份R2整體優于干旱年份，且2個年份的NDVI_降尺度TRMM數據和EVI_降尺度TRMM數據與原始TRMM數據相比，R2均有所下降，但NDVI_降尺度數據R2更接近；BIAS，RMSE，MAE均有所上升，NDVI_降尺度TRMM數據精度要好于EVI_降尺度TRMM數據。

圖6 TRMM降水值、預測降水值、殘差值及降尺度降水值

4 結 論

本文以長江流域為研究區，以2001—2019年的MODIS和TRMM產品為數據源，在充分考慮TRMM 3B43降水數據在長江流域的適用性的基礎上，基于GWR模型，選用不同的降尺度參數，對TRMM 3B43數據進行降尺度研究。并對TRMM原始數據、降尺度降水數據與氣象站點實測數據進行結果檢驗和對比分析。得出以下結論：

(1) TRMM 3B43降水數據在長江流域年時間尺度上具有較好精度，其與氣象站點實測降水數據的R2為0.911 6，且站點相關系數大多在0.75以上，表現出一定的適用性。

表3 典型年份降尺度結果對比

(2) 選取研究區5 367個樣點，采用均方根誤差(RMSE)、決定系數(R2)作為驗證指標。通過交叉驗證法對插值方式進行選擇，經分析選擇反距離權重插值法進行殘差插值結果最優。

(3) TRMM 3B43數據在兩種降尺度處理前后的空間分布特征上都整體趨于一致，表現為由東向西遞減的特點，原始TRMM降水數據局部特征較模糊，降尺度數據細節特征表現更好。NDVI_降尺度TRMM數據在與氣象站點實測數據之間的R2，BIAS，RMSE，MAE指標上都優于EVI_降尺度TRMM數據，且2006年、2010年、2012年、2015年、2019年各項指標比原始TRMM數據更好，可以反映長江流域真實的降水信息。

(4) 整體而言，兩種降尺度降水數據區間范圍相比原始TRMM降水數據區間范圍擴大，濕潤年份降尺度結果優于干旱年份。