基于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降水数据区间范围扩大，湿润年份降尺度结果优于干旱年份。