降水是反映地表环境状况和全球水循环的关键参数，是流域中水分循环和能量交换的重要组成部分[1-2]，也是表征气候变化的重要指标，同时是一切现有水文模型最关键的输入场[3-4]，高精度、高时空分布的降水观测值可有效反映地表环境状况和全球水循环的时空分布格局，有助于提高水文模型的模拟精度并据此获取更为精准的模型参数．降水观测数据一般有地面观测、雷达估测及卫星观测反演[5]三大来源．然而，现有的实测和预报降水资料均有不同的时空误差[6-9]，无法满足具体的业务需求．其中，地面观测降水的精度高但空间分辨率低；雷达降水的时空分辨率高但测量空间范围有限且易受复杂地势的影响；卫星降水产品的测量范围大且时空分辨率高，但在点上的精度不如地面观测降水．而降水预报产品的性能在各个流域的表现不一．因此，耦合具有不同时空分布特征的多源降水数据观测的同化融合技术方法成为提升降水资料精度的有效途径和当今水文研究的关注热点[10]．传统的多源降水融合方法主要基于加权平均、回归、分析等数学思想，对降水产品的误差进行处理，主要使用的方法有概率匹配方法[11]、客观统计分析法[12]、贝叶斯校正方法[13]、地理加权回归(geographically weighted regression，GWR)[14]等．以上方法通常建立在很强的假设下，且一般只考虑空间或时间因素，没有同时加入时间和空间的影响，随着计算机技术的发展和降水产品数量的增加，人工智能算法因对处理大数据有独特的优势而应用于降水数据融合和预测中[15]．其中，文献[16]提出的ConvLSTM网络既能考虑降雨序列在时间上的相关性，又能考虑降雨的空间分布特点，在短临降雨预报领域得到了广泛应用．本研究引入ConvLSTM网络进行降水融合，一方面探究ConvLSTM网络在降雨融合领域的适用性，另一方面获取长序列、高精度的实测降雨资料．基于此，提出了基于ConvLSTM网络的多源降雨融合方法．首先，获取构建模型所需的多源数据产品，包括卫星降雨、地形数据等，并对原始数据进行预处理，对卫星数据和雨量站观测数据中的缺失值和异常值进行处理，并将所有网格数据重采样至统一的空间分辨率；然后，提取网格数据构建训练样本集，并对ConvLSTM网络的参数进行训练；最后，以雅砻江流域为例，对训练得到的模型进行精度检验．1 研究区域和数据1.1　研究区域简介研究区域为雅砻江流域，雅砻江是金沙江第一大支流，位于青藏高原东部，地理位置界于北纬26°32′~33°58′，东经96°52′~102°48′之间，发源于青海省玉树县境内的巴颜喀拉山南麓，干流河道全长1 570 km，流域面积约1.3×106 km2，占金沙江(宜宾以上)集水面积的27.3%．流域南北跨越七个多纬度，且域内地形地势变化悬殊，使本流域气候气象条件在南北及垂直方向上都有明显的差异．受地形因素影响，雅砻江中上游地区降水观测网络稀疏，无法对降水进行精确观测．雅砻江流域如图1所示．10.13245/j.hust.220816.F001图1雅砻江流域示意图1.2　数据来源本研究使用GPM卫星降水数据、地面雨量站数据和数字高程模型(DEM)．根据反演算法的不同，GPM能够提供3种级别的遥感数据产品，其中三级IMERG (integrated multi-satellitE retrievals for GPM，GPM多卫星融合反演降水数据)产品是由校准后的微波所生成的红外降水估计，此外，还融合了地面观测数据，目前已更新至V06B版本．IMERG产品中的Final Run质量最高，最适合于科学研究，该产品的空间覆盖范围为60°S~60°N，空间分辨率为0.1°×0.1°，时间分辨率为日．本研究使用的数据长度为2000年6月—2020年12月，来源于NASA Earthdata (https：//earthdata.nasa.gov/)．DEM是目前用来描述流域地形地貌信息的主要手段．本研究所用的DEM数据来自于地理空间数据云(http：//www.gscloud.cn/)提供的GDEMV2 30 m分辨率原始高程数据．地面气象站点观测数据来自国家气象科学数据中心(http：//data.cma.cn/)发布的《中国地面气象日值数据集(V3.0)》，包含了中国699个基准、基本气象站1951年1月以来本站气压、气温、降水量、蒸发量、相对湿度、风向风速、日照时数和0 cm地温要素的日值数据．考虑到雅砻江流域内部分气象站点数据缺测及气象站较少的情况，本研究选用了雅砻江流域内及周边的18个气象站点与GPM数据同期的日降水观测值．各气象站点的基本信息如表1所列，空间分布如图1所示．10.13245/j.hust.220816.T001表1气象站点信息站名经度纬度拔海高度/m清水河97.08°E33.48°N4 415.40石渠98.06°E32.59°N4 200.00德格98.35°E31.48°N3 184.00甘孜100.00°E31.37°N3 393.50色达100.20°E32.17°N3 893.90道孚101.07°E30.59°N2 957.20新龙100.19°E30.56°N3 000.00理塘100.16°E30.00°N3 948.90稻城100.18°E29.03°N3 727.70康定101.58°E30.03°N2 615.70木里101.16°E27.56°N2 426.50九龙101.30°E29.00°N2 925.00越西102.31°E28.39°N1 659.50盐源101.31°E27.26°N2 545.00西昌102.16°E27.54°N1 590.90华坪101.16°E26.38°N1 230.80攀枝花101.43°E26.35°N1 190.10会理102.15°E26.39°N1 787.301.3　数据预处理本研究使用了多种不同的数据，在输入ConvLSTM网络模型训练前须要进行预处理．首先，对GPM卫星数据和雨量站观测数据中的缺失值和异常值进行处理．对于不同类型的缺失值，分在空间和时间纬度上进行线性插值．然后，为了获取特定分辨率的降水空间分布及在更小的区域内能获取更多的数据进行卷积运算，同时为了尽可能保留原值以避免引入新的误差，使用最近邻插值法将GPM数据降尺度到0.05°×0.05°．DEM数据同样也重采样到0.05°×0.05°，并根据重采样后的DEM数据用Arcgis软件提取ConvLSTM网络输入所需的高程、坡度、坡向等地表辅助变量数据．其次，考虑到不同数据具有不同的量纲，在神经网络训练时会产生影响，故对其分别进行最大值-最小值归一化，将范围限定在[0，1]．最后，提取网格数据得到训练数据．由于要考虑降水的空间信息，对于每一个地面雨量站，在卫星网格数据上以其为中心提取一个7×7(约38.5 km×38.5 km)的子网格代表当前站点的降水空间分布信息．与此相同，在整个时间序列的每个时刻提取每个气象观测站点对应的子网格，建立卫星网格数据和地面观测数据时间和空间对应的训练数据．卫星子网格提取方式示意图如图2所示，在t时刻对于第i个气象站，以气象站为中心，在GPM网格降水数据中提取一个7×7的子网格，即可获得第i个气象站周边不同时刻tk (k=0，1,⋯,n)的降水空间分布信息．10.13245/j.hust.220816.F002图2卫星子网格提取方式示意图2 ConvLSTM网络模型2.1　ConvLSTM原理介绍ConvLSTM网络是对长短时记忆网络的改进，不仅具有长短时记忆网络的时序建模能力，还能像卷积神经网络一样刻画局部特征，具备时空特性．传统的长短时记忆网络由输入门、遗忘门、状态门、输出门构成．各部分之间的关系可表示为it=σ(Wxixt+Whiht-1+Wci∘ct-1+bi)；ft=σ(Wxfxt+Whfht-1+Wcf∘ct-1+bf)；ct=ft∘ct-1+it∘tanh(Wxcxt+Whcht-1+bc)；ot=σ(Wxoxt+Whoht-1+Wco∘ct+bo)；ht=ot∘tanh(ct)，式中：it为输出门状态保留概率；σ为Sigmoid激活函数；xt为t时刻输入；∘为Hadamard乘积；ft为遗忘门状态保留概率；ct为t时刻单元状态；ot为t时刻输出门输出概率；ht为t时刻隐含层输出；W*i和bi分别为输入门的权重和阈值；W*f和bf分别为遗忘门的权重和阈值；W*c和bc分别为状态门的权重和阈值；W*o和bo分别为输出门的权重和阈值．传统长短时记忆网络内部门之间是依赖类似前馈式神经网络计算的，可以很好地处理时序数据，但是无法刻画空间数据局部特征．ConvLSTM网络将传统长短时记忆网络中输入门-状态门，状态门-状态门之间的前馈式计算换成卷积的形式，内部结构示意图如图3所示．10.13245/j.hust.220816.F003图3ConvLSTM网络内部结构示意图2.2　模型构建基于预处理后的雅砻江流域地面站点降水观测数据，GPM卫星降雨数据及经纬度、坡度、坡向、坡长等地表辅助变量数据，采用ConvLSTM网络构建多源降雨数据融合模型，Pp,t=fConv(Gt-1, Gt-2, ⋯, Gt-k, lat, lon, e, a, s)，式中：t为时间；Pp,t为t时刻指定站点融合后的降雨值；Gt-k为GPM卫星降雨值，其中k为模型须输入t时刻前k个时段的GPM降雨数据，取k=10．lat和lon分别为站点纬度和经度；e，a，s分别为气象观测站点周边网格的高程、坡向和坡度．所有输入因子均为二维矩阵．构建的ConvLSTM网络模型结构如图4所示．10.13245/j.hust.220816.F004图4ConvLSTM网络模型结构示意图2.3　精度评价为了定量评估所提出的降雨融合模型性能，使用5折交叉验证将18个雨量站数据分成5份，其中4份用于训练，剩下的1份用于测试，重复5次直到每一份都经过测试，得到5份测试结果，其中每一个雨量站都将经过测试．在整个研究区域，使用面雨量的精度指标作为评估结果．根据前人研究可知，泰森多边形由于考虑了各雨量站的权重，且当测站固定不变时，各测站的权重也不变，相对较合理，精度也较高，因此实测面雨量和融合面雨量由泰森多边形求得[17]．评估标准选择均方根误差RMSE(RRMSE)、平均绝对误差MAE(M)和相关系数(R)，计算公式为RRMSE=∑i=1n(Pp,i-Pobs,i)2/n；M=1n∑i=1nPp,i-Pobs,i；R=∑i=1n(Pp,i-P¯p)(Pobs,i-P¯obs)∑i=1n(Pp,i-P¯p)2(Pobs,i-P¯obs)2，式中：n为降水数据序列长度，Pp,i和Pobs,i分别为第i个时段融合降雨面雨量和实测面雨量；P¯p和P¯obs分别为融合降雨面雨量和实测面雨量的均值．2.4　建模步骤基于ConvLSTM网络进行GPM卫星降水数据和气象站点数据的融合，建模的主要步骤如下．步骤1 收集并整理雅砻江流域地面站点降水观测数据，包括GPM卫星降雨数据、站点经纬度，坡度、坡向、坡长等数据，构建输入样本集．步骤2 对输入样本集进行归一化处理．步骤3 根据输入样本集构建ConvLSTM网络模型，采用5折交叉验证方法划分训练集样本和测试集样本，对网络参数进行训练．步骤4 通过泰森多边形将点雨量转换为面雨量，利用均方根误差、平均绝对误差和相关系数对模拟结果进行精度评定．3 结果分析3.1　融合降雨精度评估为分析ConvLSTM网络降雨融合精度模型，分别使用ConvLSTM网络模型、GWR方法和LSTM模型进行降雨数据融合，并将计算结果与GPM卫星原始数据精度进行对比，表2为不同降雨融合模型精度．图5为不同降雨融合模型逐日面雨量(Pd)与观测站点面雨量(Ps)相关性．其中，由GWR方法、LSTM模型、ConvLSTM网络模型和气象观测站点得到的雨量值采用泰森多边形法转换为雅砻江流域逐日面雨量值，GPM数据采用流域内所有网格点的雨量均值作为面雨量值．10.13245/j.hust.220816.T002表2不同降雨融合模型精度模型RRMSE/mmM/mmR原始GPM5.151.350.735GWR4.651.270.755LSTM1.970.820.911ConvLSTM1.800.790.91110.13245/j.hust.220816.F005图5不同降雨融合模型逐日面雨量与观测面雨量相关性由表2可知：LSTM模型和ConvLSTM网络模型融合降雨结果精度较原始GPM卫星降雨数据和GWR方法融合降雨结果精度均有明显提升．其中，ConvLSTM网络模型融合降雨结果精度提升更为明显．原因在于，ConvLSTM网络模型考虑了输入因子在空间上的分布特性，融合了气象观测站点周边的地形信息，而LSTM网络须要将在空间上呈二维分布的输入因子展开成一维向量进行输入和计算，失去了输入因子的空间分布特征．此外，GWR未考虑降雨在时间上的相关性，输入因子与输出因子的时间为同一时刻，且仅考虑了气象观测站点与研究区域内网格点的距离因素，在计算中将一些距离观测站点较远的网格点信息也输入到模型中，引入了一些非必要的输入信息，从而影响计算结果的精度．综上所述，ConvLSTM网络模型能够同时考虑输入因子的空间分布特性和时间上的相关性．为进一步分析不同融合降雨与站点观测降雨趋势，分别依据GWR方法融合降雨结果、ConvLSTM网络模型融合降雨结果和气象观测站点雨量值，采用泰森多边形法计算了雅砻江流域逐月面雨量值，GPM数据则采用流域内所有网格点的雨量均值作为面雨量值．不同降雨融合模型逐月面雨量与观测面雨量对比如图6所示，图中P为逐月面雨量．由图6可知：四种降雨数据的变化趋势基本符合，但GPM数据和GWR方法融合降雨结果在峰值处较由气象观测站点雨量求得的面雨量值偏高，而LSTM模型、ConvLSTM网络模型融合降雨结果在峰值处则与气象观测站点面雨量更为一致且ConvLSTM网络模型对月尺度面雨量峰值的模拟效果更好．由GPM数据求得面雨量峰值偏高的现象与其他学者一致[18]，说明本研究提出的ConvLSTM网络模型能够在一定程度上改善GPM数据对于雅砻江面雨量峰值的估计．10.13245/j.hust.220816.F006图6不同降雨融合模型逐月面雨量与观测面雨量对比3.2　融合降雨空间分布特征为分析模型对降雨空间分布的模拟能力，研究采用ConvLSTM模型进一步构建雅砻江流域0.05°分辨率的日降雨融合数据集．研究选取了2018年6月13日的降雨模拟结果进行分析，该日雅砻江流域实测面雨量较大(20.26 mm)，对流域防洪造成较大压力．依据所选取的18个气象站点的同期雨量，采用反距离插值法绘制实测降雨的雨量分布图，并将ConvLSTM模型融合降雨的分布与实测降雨分布进行对比分析．ConvLSTM模型融合降雨和实测降雨分布如图7所示．10.13245/j.hust.220816.F007图7融合降雨和实测降雨分布由图7可知：融合降雨结果能够捕捉到雅砻江下游西南角的暴雨中心，但对于暴雨中心雨量等级存在明显的低估，主要是由于作为关键输入数据的GPM卫星数据本身对于极端天气事件的监测不足，对于强降水存在明显低估．另外，ConvLSTM降雨融合数据可以把握全流域的降雨总体分布情况，且相对原始GPM数据拥有更高的空间分辨率，可为流域水文研究提供相对可靠的数据来源．4 结论本研究以雅砻江流域为研究区域，以地面站点观测降雨为因变量，结合GPM卫星降雨、经纬度、高程、坡度、坡向为输入因子，基于ConvLSTM网络构建了多源降雨数据融合模型，并将所提模型的计算精度与原始GPM卫星数据、GWR方法和LSTM模型计算精度进行了对比．进一步应用所提模型构建了雅砻江流域0.05°分辨率的日降雨融合数据集，并以2018年6月13日典型降雨为例，分析了多源降雨数据融合模型结果的空间分布特征，研究结论如下．a. 基于ConvLSTM网络构建的多源降雨数据融合模型能够在保证降雨量精度的同时，提升降雨数据的空间分辨率，且ConvLSTM网络模型的融合降雨精度高于原始GPM卫星数据、GWR方法和LSTM模型计算精度，能够一定程度上改善GPM卫星数据对于降雨峰值的估计．b. 基于ConvLSTM网络构建的多源降雨数据融合模型充分考虑了地面降雨与卫星降雨及地形特征的非线性关系，能够较为准确地模拟出流域降雨的空间分布情况，且能够准确展现暴雨中心位置，但对于日尺度暴雨中心雨量等级的估计仍有一定的提升空间．