基于集合卡尔曼滤波法的清江流域多源融合降水分析

郭家力，郭东淏，丁光旭，李 颖1，，李英海1，，杨 旭，张海荣

（1.水电工程施工与管理湖北省重点实验室（三峡大学），湖北 宜昌 443002； 2.三峡大学水利与环境学院，湖北 宜昌 443002；3.智慧长江与水电科学湖北省重点实验室（中国长江电力股份有限公司），湖北 宜昌 443133）

0 引 言

降水是最基本的气象要素之一，准确掌握流域降水的时空分布特征对实现数字孪生流域“四预”功能具有重要意义［1-3］。目前获取降水资料的主要途径主要有3种：地面雨量站、地基天气雷达估测和气象卫星反演［4］，各自优缺点十分明显，但在一定程度上可以互补［5］。地面雨量站对降水量的观测是获得降水数据最直接也是最常用的方法，在整个降水观测体系中起到基础性的作用，可以获得高精度的单点数据，但地面站点观测无法覆盖到无人区以及地形相对复杂的区域，从而限制了站点观测数据的使用；气象卫星反演降水通过间接手段进行观测，不受地理条件的限制，具有全天候、全覆盖的特点，然而其精度会受到传感器性能、反演算法以及云层性质的影响，虽然可一定程度上反映降水的时空分布，但在单点上精度低于地面观测数据［4，6］；雷达通过间接手段观测降水，其优点在于能够跟踪雨区范围，在中小尺度地区的短期强降水监测中具有较强的优势，但雷达对降水的观测范围有限，且易受地形等因素影响，估测的降水数据误差较大［7］。受制于降水类型、地形地貌、反演算法等因素，单一降水产品难以准确反映区域降水过程，导致各降水产品在不同区域有着不同的性能。

因此，如何通过多源降水数据融合的方式研制高质量和高时空分辨率的降水产品，一直是水文气象领域研究的热点问题之一［8］。王筱译［9］等将加权最小二乘估计法应用于降水数据融合，引入SM2RAIN-CCI降水数据对IMERG数据进行校正，并以渭河流域2014年至2015年的地面实测降水数据进行效果检验，融合数据的各项统计误差较原始IMERG数据均有较大的改善；谭伟伟［10］等提出了点面融合方法和站点偏差校正估计两个多源数据融合方案，考虑了经纬度、DEM、IMERG插值数据等辅助变量，融合了湖北省2016年7月19日的IMERG日尺度降水数据与气象站点数据资料，有效提升了IMERG降水数据的精度；高真［11］等以地面雨量站降水数据作为参考，评估了NLDAS2、Stage IV、TMPARP、TMPART和IMERG降水数据的精度，并分别驱动分布式水文模型，发现多源遥感降水数据融合产品精度更高，更适合于水文模拟应用；孟庆博［9］等基于地面观测数据和TRMM、CMORPH、CHIRPS、PERSIAAN-CDR、GL‐DAS_Noah多源卫星降水数据，使用了集合卡尔曼滤波算法进行数据融合同化，同化后的5种降水产品精度均有显著提升，证明集合卡尔曼滤波算法在降水数据融合同化中有较大潜力。从上述研究可以看出，目前在多源降水融合研究方面，集合卡尔曼滤波被证明是一种具有较大潜力的方法，但目前的应用研究却开展极少。

基于此，本文将基于集合卡尔曼滤波算法融合多源降水资料，生成一套高精度的降水融合产品。首先，收集建立融合模型所需的多源数据，包括地面站点地形数据、实测数据和卫星降水数据，然后对原始数据进行预处理，将所有网格降水数据重采样至统一的空间分辨率，处理成融合模型所需要的数据，最后使用集合卡尔曼滤波算法对多源降水数据进行融合，得到一套0.05°×0.05°的融合降水数据。

1 研究区域概况与数据获取

1.1 流域概况

清江属于长江的一级支流，位于东经108°58'～111°21'，北纬29°41'～30°87'之间，属长江上游地区，流域位于湖北省境内。清江干流发源于恩施州利川市境内，全长423 km，平均比降3.2%，流域面积1.67 万km2。清江流域内部群山耸立、河流深切、高差悬殊，高山地区（海拔>1 200 m）占流域总面积的28%；中高山（海拔800～1 200 m）面积为1.2万km2，占52%；低山丘陵区（<800 m）为4 652 km2，占20%，流域最大相对高差2 217 m。清江流域雨量丰沛，水汽主要来自太平洋季风和印度洋季风，年平均降水量为1 561 mm，最大年降雨量达1 700多 mm，最少年降水量仅900 mm左右，其中80%以上的年降水量集中在汛期5-10月，各月降水差异明显，属于亚热带季风气候。清江流域地理位置如图1所示。

1.2 数据来源

1.1.1 地形高程数据

地形高程数据选择SRTM DEM 90 m分辨率原始高程数据，数据从地理空间数据云下载（http：//www.gscloud.cn/）。

1.1.2 地面站点数据

本文收集到了流域内部28个及流域周边5个共33个地面观测站点的逐日降水数据，站点分布如图1所示，所选数据通过了严格的质量控制，资料完整，数据可靠，最终得到1998年1月1日-1999年12月31日的逐日降水资料。在日尺度进行多源降水数据融合，数据序列长度730 d，基本满足融合结果的精度检验要求。

1.1.3 卫星降水数据

本文所用的卫星降水数据包括TRMM（Tropical Rainfall Measuring Mission）和CMORPH（Climate Prediction Center morph‐ing technique）数据，其空间分辨率均为0.25°，时间分辨率为1 d。TRMM卫星由美国国家航空航天局和日本宇宙航空研究开发机构联合研制，是专门用于定量测量热带、亚热带降雨的气象卫星，现已扩展到近似覆盖全球范围［13］，研究采用美国国家航空航天局发布的TMPA（TRMM Multi-satellite Precipitation Analysis）第7版日降水数据。CMORPH数据是美国国家海洋和大气管理局发布的卫星降水数据集［14］，本文采用的是CMORPH_V1.0_BLD日降水数据。

1.1.4 再分析降水数据

为了进一步验证融合数据集的可靠性，选用CMFD（China Meteorological Forcing Dataset，CMFD）、ERA5（ECMWF Reanaly‐sis v5）和MSWEP（Multi-Source Weighted-Ensemble Precipita‐tion）三种再分析数据集与其进行对比。CMFD是中国科学院青藏高原研究所开发的近地面气象与环境要素再分析数据，该数据以国际上已有的GLDAS（Global Land Data Assimilation Sys‐tem）等再分析数据和TRMM数据为背景数据，融合中国气象局常规气象观测资料制作而成［15］。ERA5是欧洲天气预报中心ECMWF（European Centre for Medium-Range Weather Forecasts）对全球气候的第五代大气再分析数据集，它使用了ECMWF综合预报系统中的四维变分同化系统和模式预报系统，降水部分融合了TRMM等数据集、雷达降水数据集和地面气象观测资料，时间分辨率为6 h，空间分辨率为0.25°×0.25°［16］。MSWEP是美国普林斯顿大学研制的融合降水数据集，该产品基于Budyko框架和全球径流观测对地形复杂区域降水进行校正，较大程度上减小了估算误差，并进一步融合了CMORPH等卫星降水数据、ERA-Interim（ECMWF Reanalysis Interim）等再分析数据和地面气象观测资料［17］。以上多源降水数据总结于表1。

表1 地面站点、卫星和再分析降水数据基本信息Tab.1 Basic information on ground stations，satellites and reanalysis of precipitation data

2 研究方法

2.1 空间插值方法

如表1所示，考虑到地面站点和卫星降水数据的空间分辨率不一致，在数据融合之前需要在网格上统一。首先通过薄板光滑样条模型（ANUSPLIN）插值方法将分布不规则的站点观测数据插值为0.05°×0.05°分辨率的规则网格数据。ANUSPLIN是一种使用薄盘光滑样条函数进行空间插值和分析的工具，在气象变量插值中取得了优异的效果，广泛应用于气象数据集的构建。除了独立的样条变量外，ANUSPLIN允许引入参数化的线性子模型或线性协变量子模型［18］。

清江流域地形复杂，海拔差异大，而地形对于降水的动力、热力和微物理过程具有重要的影响，是导致区域降水量变化的一个主要因素［19］。因此选取经度和纬度作为独立变量，海拔作为协变量，对站点观测数据进行插值，得到的结果将其规定为插值数据SI。

2.2 融合方法

集合卡尔曼滤波算法是集合预报与卡尔曼滤波方法的结合，它是贝叶斯信息估计理论的一种直接实现［20］。集合卡尔曼滤波算法通过蒙特卡罗方法计算状态的预报误差协方差，用集合的思想解决实际应用中背景误差协方差矩阵估计和预报困难的问题，有效降低了计算量。与标准的卡尔曼滤波相似，集合卡尔曼滤波算法主要分为预测和更新两个步骤［21］。

在k=0时刻对状态变量集合进行初始化（本文仅有降水这一个状态变量），将k时刻的状态变量代入模型，得到状态量k+1时刻的预测值。

计算卡尔曼增益K，并对k+1时刻的观测值对集合的状态进行更新得到分析值。

将地面站点插值数据作为融合算法观测值，将卫星降水数据作为融合算法背景值，统一在0.05°×0.05°分辨率下的网格进行数据融合。

2.3 精度评定方法和指标

为了评估融合模型的精度，并对比主流的再分析数据，采用留一交叉验证法进行验证。留一交叉验证（Leave-One-Out Cross-Validation，LOO-CV）是评估模型泛化能力的方法，具体做法为：每次使用5个流域周边气象站点和清江流域内28个气象站点中的27个站点作为模型输入，剩下的1个气象站点留做精度评定，如此逐一循环28次，使得流域内每个气象站都参与了精度评定，得到插值数据在流域内每个站点的交叉验证值。卫星数据、融合数据以及再分析数据则直接提取28个站点所在网格的数据进行精度评定。使用相关性系数R，平均绝对误差MAE和均方根误差RMSE作为具体的评价指标，其表达式如公式（8）～（10）所示。

式中：Gi表示第i天的站点实测降水；Si表示第i天的卫星降水（或插值降水、融合降水和再分析降水），mm；N表示序列长度，d。

3 结果与分析

3.1 融合模型精度评价

通过采用留一交叉验证的方法，得出清江流域内28个站点日尺度下各评价指标。图2为背景数据（TRMM、CMORPH和SI）、再分析数据（CMFD、ERA5和MSWEP）与MSAP融合数据评价指标箱线图。其中，红色代表3种背景数据，绿色代表MSAP数据，蓝色代表3种再分析数据。

图2 日尺度融合数据相关性系数、平均绝对误差、均方根误差的箱线图Fig.2 Boxplots of R， MAE， and RMSE of daily-scale fusion data

首先，对比背景数据和MSAP数据的精度。图2（a）为相关性系数箱线图：TRMM、CMORPH和SI的中位数分别为0.67、0.81和0.88，MSAP中位数为0.89，相较于TRMM、CMORPH和SI分别提高了32.1%、9.3%和1.1%；数值分布上，MSAP在各站点均大于0.8，而SI有3个站点小于0.8，可见MSAP略优于SI。图2（b）为平均绝对误差箱线图：TRMM中位数为3.7 mm/d，最小值为3.1 mm/d；CMORPH中位数为2.5 mm/d，数据分布于1.4～3.4 mm/d之间；SI中位数为1.9 mm/d，最大值为3.6 mm/d；MSAP中位数为1.9 mm/d，与SI持平，但优于TRMM和CMORPH；数据分布上，MSAP最大值为2.9 mm/d，低于SI，数据分布均匀，可见MSAP在平均绝对误差指标上能基本保持与三种背景数据中的最优者SI持平。图2（c）为均方根误差箱线图：TRMM数据均在7 mm/d以上，误差较大；CMORPH分布于4.0～10.9 mm/d之间，中位数为6.9 mm/d；SI分布于2.8～12.74 mm/d之间，中位数为5.5 mm/d，多个站点值大于8 mm/d；MSAP相较于TRMM和CMORPH分别提高了41.2%、24.3%，与SI持平；MSAP的均方根误差数值分布上极差最小，集中分布于3.7～9.5 mm/d之间，最大值为9.5 mm/d，均小于3种背景数据的最大值。通过对3个评价指标对比结果的定量分析可以得出结论，EnKF融合算法融合了各背景数据的优势，MSAP数据的精度优于3种背景数据。

其次，对比MSAP与现有再分析降水数据的精度。从图2绿色、蓝色箱线图对比可以看出，无论是相关性系数、平均绝对误差还是均方根误差，MSAP与3种再分析降水数据中的CMFD最为接近。MSAP相关性系数中位数为0.89，而CMFD为0.78，整体相关性系数优于CMFD；尽管CMFD的平均绝对误差的中位数小于MASP，但其分布范围更宽，说明CMFD空间异质性更大，而均方根误差方面则是MASP明显优于CMFD。MSWEP的分布在3个指标上略逊于CMFD，且无论是中位数还是分布范围二者均类似。ERA5在3个指标中表现均为最差，精度远低于CMFD、MSWEP和MSAP。总体上，MSAP在各个指标上全面优于3种再分析数据，精度排序为MSAP > CMFD > MSWEP >ERA5，进一步说明了EnKF算法在降水数据融合方面的有效性。

3.2 降水精度评价指标空间分布

箱线图仅能表达降水精度评价指标误差的范围，但并不能直观展示评价指标的空间分布状况。为此，图3展示了背景数据、再分析数据和融合数据的精度评价指标空间分布图。由图3背景数据（a）～（c）和融合数据图3（g）可知，在3种背景数据各有所长。SI数据整体精度高，但空间异质性较大，特别是在咸丰、汪营、五峰等流域边界地区精度较低，这是因为融合产品质量与站网密度呈正相关关系［22］；CMORPH各站点精度空间分布较为均匀，但在流域上游地区的精度低于其他地区，且在当阳坝地区均方根误差偏大；TRMM在总体精度和空间分布上均不如前面两种数据，但从平均误差指标来看，流域内的部分站点如西流水、茅田等精度高于SI和CMORPH。MSAP综合了卫星数据和高精度插值数据，改善了背景数据流域边界地区部分区域精度低的缺点：在当阳坝等站点，MSAP精度高于卫星数据，与SI精度持平；在五峰和咸丰等站点，MSAP精度高于SI，略低于CMORPH；在汪营等站点，MSAP精度则优于3种背景数据，证明了EnKF算法的有效性。在图3（d）～（f）的再分析数据中，CMFD与MSWEP各站点精度评价指标的空间异质性高于MSAP数据，两种数据在流域下游的精度高于流域上游。CMFD在所有站点的精度评价指标都低于MSAP，MSAP融合了清江流域自动站观测资料，而CMFD融合了中国气象局常规气象观测数据，MSAP的背景数据在站网密度上高于CMFD，这是MSAP精度高于CMFD的重要原因［22］。

3.3 场次强降水量空间分布分析

基于更加直观地验证本文建立的融合模型所得到的降水与3种再分析降水产品的空间分布特征，根据融合后的降水对清江流域1998-1999年间的不同历时的场次强降水进行分析，以适应不同场景下降水产品的应用需求。

1998年长江上游从6月中旬开始到8月中旬长达两个多月的强对流云团的频繁发生是1998年长江上游洪峰多，中游高水位长时间维持的重要原因［23］。这一时期长江上游地区暴雨洪水频发，清江流域自5月中旬以来，也连续发生了8次洪峰，其中第四、五、七、八次洪峰与长江干流第一、二、四、六次洪峰不同程度地遭遇［24］。选取四次遭遇中洪峰最大的第四次和第八次进行降水空间分析，其历时分别为5日和2日，对应的流域平均面雨量为186和86 mm。

1998年6月28日-1998年7月2日清江流域发生了1998年最大的一次强降雨，此次降雨持续5天，暴雨中心位于流域北部的茅田地区，单站降水量达到330 mm，流域平均面雨量达到了186 mm，此次强降水也引发了清江流域1998年最大的一次洪水（洪峰流量12 300 m3/s，重现期20年一遇）［24］。图4给出了各数据集对此次降雨的捕捉情况。可以看到各数据集都捕捉到了这个全流域强降水的特性，5日累计降水量都在100 mm以上，但在降水空间分布上有几点差异：CMFD、MSWEP和MSAP能反映出暴雨中心，但是显示暴雨中心位于流域北部，并向南部递减，在流域南部上游地区达到最低；ERA5趋势与CMFD趋势一致，降水量最小；MSWEP与其他数据趋势一致，降水量级为最大；MSAP显示暴雨中心位于流域西北部，并向东南递减，在量级上与CMFD保持一致。

清江流域属于亚热带季风气候，夏季易受到东南季风的影响。1998年8月15日-1998年8月16日清江流域出现了一次强降雨，暴雨中心位于流域东南地区的五峰站，48 h降水量为146 mm，流域平均面雨量为86 mm。图5给出了各数据集对此次降雨的捕捉情况。从图5中可以看出，各数据集均观测到了此次降雨的空间分布：ERA5与其他3种数据有较大差异，其降雨中心位于流域中上游；CMFD、MSWEP和MSAP的降雨中心均位于流域东南部，并向流域上游递减，其中CMFD与MSAP的空间分布更相似。

图5 1998年8月15日-8月16日典型日降水分析Fig.5 Analysis of typical daily precipitation from August 15 to August 16， 1998

4 结 论

在日尺度上通过集合卡尔曼滤波法（EnKF）融合算法对清江流域1998～1999年的地面插值数据及2种卫星降水产品进行数据融合，得到了清江流域的0.05°×0.05°融合降水序列，利用留一交叉验证的方式，验证了融合产品的精度。同时采用了国际上主流的再分析数据与本文融合降水数据进行对比，并对清江流域内的场次强降水空间分布进行了分析。

（1）交叉验证结果显示，EnKF算法改善了卫星降水整体精度低和插值数据流域边界精度低的缺点，使得MSAP数据普遍优于融合前的单一降水产品；融合数据同样优于3种再分析数据，CMFD数据与融合数据最为接近，MSWEP数据次之，ERA5数据精度最差。

（2）融合产品能准确捕捉清江流域场次强降水的空间分布特征，MSAP融合数据与我国主流融合降水CMFD数据呈现相似的空间分布，且强降雨历时越长，空间分布特征越接近。