引言

降水是陆面、水文和生态模型中的重要输入参数之一，其时空分布和变化对水文循环、地-气交互和水资源利用等具有重要影响。近年来，随着观测技术迅速发展，气象观测数据已从单一地面、高空观测发展到一体化地基、空基、天基等多源观测系统，进而有助于获取时空分辨率更高、覆盖范围更广的多源降水观测资料。然而，任何一种观测手段都存在优势和缺陷，如地面自动站、雷达探测可以获取准确的降水资料，但存在不能实现全覆盖等问题(Duan et al., 2016)；卫星降水资料具有覆盖范围广、全天候连续观测的优势，但并非直接实际观测，产品精度也较低(刘元波等，2011)。综合应用多种平台观测资料优势，融合多种来源的地面、卫星、雷达等各类降水观测资料是解决单一降水资料不确定性等问题的有效途径。不断提高资料的数据精度已成为高质量降水产品研制的发展趋势。目前，国际上已发展了较多基于台站观测、卫星反演和雷达估测的高分辨率融合降水产品，如GPCP (Global Precipitation Climatology Project)月产品(Huffman et al., 1997)、GPCP旬产品(Xie et al., 2003)、GPCP 1DD (One-Degree Daily)日产品(Huffman et al., 2001)、TRMM 3B42 (Tropical Rainfall Measuring Mission 3B42) (Huffman et al., 2007)降水产品等。国内气象学者也开展了较多的多源降水融合产品研制工作。Shen等(2014)利用“概率密度匹配+最优插值(PDF+OI)”两步融合方法，对中国地面站点观测降水与卫星降水产品(FY、CMORPH)进行融合，制作中国区域10 km分辨率逐小时融合降水产品。潘旸等(2015)采用“概率密度匹配+贝叶斯模型平均法+最优插值(PDF+BMA+OI)”方法，对中国地面站点观测降水、卫星降水和雷达估测降水进行融合分析，制作中国区域5 km和1 km分辨率逐小时融合降水产品。师春香等(2013)采用多重网格变分分析方法(STMAS)，制作中国气象局陆面数据同化系统(CLDAS)的大气驱动场数据(包括2 m气温、气压、湿度、风速、降水和太阳辐射等)，获得了准确的土壤湿度模拟结果。

当前，卫星降水虽可弥补站点观测降水的空间覆盖不足，但其结构和数据精度偏差较大，与之相比，雷达估测降水(QPE)在降水结构和数据精度上均有明显提高，且在区域尺度上可实现时空全覆盖，因此雷达观测资料是研制区域尺度上高质量、高时空分辨率降水产品不可或缺的数据源。随着IT技术不断发展和我国多普勒天气雷达网布设趋于完善，国内众多机构相继开展了全国或区域组网的多种雷达QPE产品研制工作。中国气象局气象探测中心可实时生成全国单站和组网的1 km·h^-1分辨率的雷达QPE产品(谷军霞等，2018)，江西省气象局可实时制作江西省组网的1 km·h^-1分辨率的雷达QPE产品(Wu et al., 2018)。上述雷达QPE产品为融合降水分析提供了重要数据源。多源数据融合方法通常依靠地面站点观测资料校正雷达、卫星等遥感估测结果(雷晓辉等，2018)，因此站点数据量及数据质量对融合效果至关重要。依托江西省山洪地质灾害防治气象保障工程及江西省山洪地质灾害防治县级非工程措施建设等项目支持，2015年江西省气象部门与水利部门开展了所属自动气象站、自动水位雨量站(以下简称“水文站”)每小时雨量观测数据共享工作，实现2 545个气象自动站和3 683个水文站信息共享，使全省雨量站密度由9 km²提升到5 km²；此外，气象站和水文站的站点布设侧重点不同，水文站主要分布在重点中小河流重要河段、山洪沟、存在病险隐患水库和大中型水闸、鄱阳湖区重点圩垸，其降水数据能有效弥补气象资料探测空间的不足。李显风等(2017)对比分析水文和气象降水数据表明，水文站与气象站观测降水的一致性较好，这是开展气象和水文雨量信息融合的前提和基础。

目前，现有降水融合产品仅融合了气象站降水资料。从融合方法上，通过融合更多站点资料可进一步提高产品的数据质量，但由于气象站和水文站所使用的采集设备及其采集原理不同，造成两套观测降水存在一定偏差，进一步融合水文站降水资料究竟对融合分析的改进效果如何，还有待验证。因此，本文利用2019年主汛期5—6月强降水时段江西省3 476个水文站、6 534个气象站(含周边省份部分站点)小时降水数据和江西省气象局研制的小时雷达估测降水数据，在对各类降水数据进行严格质控处理的基础上，采用多重网格变分分析方法，开展了1 km×1 km逐小时江西省气象站、水文站降水数据与雷达背景场降水的融合试验，并对试验结果进行了分析评估。

1 数据预处理与融合方法 1.1 数据及预处理 1.1.1 气象站观测降水

本研究融合分析区域覆盖整个江西省(113°—119°E，24°—31°N)，气象站观测降水包括江西省国家级气象观测站(以下简称国家站)、区域气象自动站(以下简称区域站)及周边6省部分国家站和区域站小时观测降水(站点分布见图 1)，共计6 534个站点，其中江西省国家站点93个。周边省份站点数据通过省际气象共享业务实时收集，所有气象站点观测降水均通过省级气象资料业务系统(MDOS)进行严格质控处理。

1.1.2 水文站观测降水

目前，国内气象站降水观测主要采用上海气象仪器厂生产的SL3-1型双翻斗雨量计，其采集分辨力为0.1 mm；水文站主要采用南京水利水文自动化研究所防汛设备厂生产的JDZ05-1型雨量计，其采集分辨力为0.5 mm (表 1)。气象、水文部门通过数据共享，双方能获得的测站降水观测资料平均增加一倍。有关研究结果(李显风等，2017)表明，水文站与气象站观测降水具有很好的一致性。但因两个部门的采集设备、采集方式不同，两套融合降水存在一定偏差，通过两者的融合应用，能否获得更高质量的实况降水是一个需要研究的问题。

对水文站观测降水资料，首先，参考全国自动气象站小时降水资料的质量控制方法(任芝花等，2010；张强等，2016)，对收集的水文站降水数据进行气候极值、区域界限值、时间一致性、空间一致性和数据完整性等快速质控处理；然后，通过与邻近气象站降水、雷达估测降水进行交互检验分析，剔除质控异常站点及异常时次数据，得到质控后3 476个全省均匀分布的水文站点(图 1)逐小时观测降水数据。

1.1.3 雷达估测降水

雷达估测降水(QPE)是区域尺度上降水融合分析的重要数据源。本研究使用的雷达估测降水为江西省气象局实时业务产品，该产品采用基于雷达回波顶高的动态反射率因子与降水关系(Z-R关系)方法进行雷达降水反演。首先，对单站雷达组合反射率(CR)和回波顶高度(ET)产品进行组网拼接；然后，利用拼接的逐时次产品进行小时累计计算；最后，根据不同ET等级(如1、2、3、4、5、6、7、8、9、10、11、12、13、14、15 km)构建每个等级的Z-R关系，利用该关系进行雷达估测降水实时计算。该方法较传统固定Z-R关系和基于反射率分级的动态Z-R关系方法对产品精度有显著改善(Wu et al., 2018)。该产品空间分辨率为1 km，时间分辨率为1 h，其空间覆盖整个江西区域(112°—120°E，23.3°—31.3°N)。

对QPE产品的质量控制，首先，采用中值滤波方法消除QPE产品因地物等杂波或雷达探测故障造成的零星小降水及异常降水，滤波器的窗口设置为3×3个网格。然后，采用概率密度匹配方法(PDF) (宇婧婧等，2013；曹萍萍等，2018)订正QPE产品的非独立性系统性偏差。

为了验证上述订正方法的适用性及订正效果，采用地面观测降水对雷达估测降水的订正结果进行检验分析。图 2给出2019年6月6日15时赣东北区域地面观测降水实况以及订正前后的雷达估测降水分布。地面观测降水显示(图 2a)，该时次有两个明显的强降水中心，其中心区雨强超过20 mm·h^-1，部分站点高达30 mm·h^-1以上。雷达降水订正前(图 2b)，其降水具有与地面观测降水相似的空间分布形态，但降水中心强度明显偏弱；订正后(图 2c)，雷达降水的空间分布形态没有改变，但降水量值与地面观测更接近。

为了进一步检验上述订正效果，图 3给出2019年6月1—10日江西省逐小时雷达降水订正前后的误差时间序列。从中看到，采用滤波和概率密度匹配方法对雷达降水进行订正后，绝大多数时刻其系统偏差减少幅度均较明显，更加接近0线，说明该方法能从一定程度上消除资料的系统偏差。

1.2 融合方法

多源数据融合是发挥日益丰富的观测数据优势的有效途径，通过对包含观测误差的多种观测数据进行综合分析以获取最优数据估计。为得到理想的降水融合效果，本文采用Xie等(2011)提出的多重网格变分分析方法(the Space and Time Multiscale Analysis System，STMAS)对站点观测数据和雷达降水背景场进行融合分析，解决最优插值和二维变分等传统融合方法对背景场误差协方差矩阵的依赖，获得更为准确的融合分析结果。

在STMAS方法中，每一重网格上的目标泛函数采用如下形式(Xie et al., 2011)

$ \begin{aligned} J^{(n)}=& \frac{1}{2} X^{(n) T} X^{(n)}+\frac{1}{2}\left(H^{(n)} X^{(n)}-\right.\\ &\left.Y^{(n)}\right)^{T} O^{(n)-1}\left(H^{(n)} X^{(n)}-Y^{(n)}\right) \end{aligned} $

(1)

式(1)中，Y=Y_obs-HX_b，X=X_a-X_b，其中，Y为观测场与背景场的差(其计算见公式2)，X代表相对模式场矢量的修正矢量(由变分数据同化系统计算得到)，H为模式网格到观测点的双线性插值算子，Y_obs为观测场矢量，X_b为背景场矢量，X_a为分析场矢量；O为观测误差协方差矩阵；n表示第n重网格；N为网格总重数。上标T和-1分别表示矩阵的转置矩阵和逆矩阵。Y⁽¹⁾、Y⁽ⁿ⁾分别表示为

$ \left\{\begin{array}{ll} Y^{(1)}=T^{\mathrm{obs}}-H^{(1)} X^{\mathrm{b}} & (n=1) \\ Y^{(n)}=Y^{(n-1)}-H^{(n-1)} X^{(n-1)} & (n=2, 3, \cdots, N) \end{array}\right. $

(2)

多重网格分析方法先从最粗的网格开始分析(n=1)，将雷达降水背景场映射到最粗的网格点上，得到最粗网格上的观测增量场X⁽¹⁾，依次类推，分别得到n=2, …, N的观测增量，进而得到最终的分析增量。最终分析结果为各重网格分析结果的叠加，即X^(a)

$ X^{\mathrm{a}}=X^{\mathrm{b}}+\Delta X=X^{\mathrm{b}}+\sum _{n=1}^{N} X^{(n)} $

(3)

根据融合输入资料和输出产品特征，本研究修改了STMAS算法的资料读写接口，实现本地资料接入与融合处理。根据研究区天气、气候与地理环境特征，采用气象站降水资料开展融合算法关键参数(如影响半径等)的融合试验。以影响半径(r)为例，分别选取3、5、10、15、20等不同影响半径(单位为网格数)，得到不同的融合分析结果，并对分析结果进行了检验评估(表 2)。表 2显示，当r达到10及以上时，相关系数(COR)、均方根误差(RMSE)、平均误差(ME)、平均绝对误差(MAE)等各项统计指标最优，随着r增大，各项统计指标差别很小，考虑到降水的空间异质性，本文r取值为10。通过开展关键参数的融合试验分析，确定了最佳参数化方案，并以改进的STMAS方法和参数化方案，开展降水融合试验研究。

1.3 融合试验

为充分检验融合试验结果，科学评价水文站降水资料对融合结果的影响，选取江西省全部国家站和遴选的部分区域站(共计339个站点)作为检验站点，因检验站点不参与融合分析，该检验方式为独立性检验。为验证进一步融合水文站对降水融合的改进效果，设计了“气象站+雷达”、“气象站+水文站+雷达”两组对比试验方案(表 3)。具体试验设计为:方案Ⅰ，采用中值滤波、PDF进行雷达估测降水订正处理，采用STMAS方法进行研究区气象站降水雷达降水数据融合分析；方案Ⅱ，采用快速质控和多源交互验证方法，对水文站降水进行质控处理，采用中值滤波、PDF进行雷达估测降水订正处理，采用STMAS方法进行研究区气象站降水、水文站降水与雷达降水数据融合分析。

1.4 评价方法

采用常用统计指标，对融合试验结果进行定量评价，统计指标包括相关系数(COR)，均方根误差(RMSE)、平均误差(ME)、平均绝对误差(MAE)、降水评分(TS)等，其统计指标计算式如下

$ C O R=\frac{\sum _{j=1}^{N}\left(G_{i}-\bar{G}\right)\left(O_{i}-\bar{O}\right)}{\sqrt{\sum _{j=1}^{N}\left(G_{i}-\bar{G}\right)^{2}} \sqrt{\sum _{i=1}^{N}\left(O_{i}-\bar{O}\right)^{2}}} $

(4)

$ R M S E=\sqrt{\frac{1}{N} \sqrt{\sum _{j=1}^{N}\left(G_{i}-O_{i}\right)^{2}}} $

(5)

$ M E=\frac{1}{N} \sum _{i=1}^{N}\left(G_{i}-O_{i}\right) $

(6)

$ M A E=\frac{1}{N} \sum _{i=1}^{N}\left(\left|G_{i}-O_{i}\right|\right) $

(7)

$ T S_{k}=\frac{\sum _{i=1}^{N}\left[G_{i} \geqslant L_{k}\right]\left[O_{i} \geqslant L_{k}\right]}{\sum _{i=1}^{N}\left(\left[G_{i} \geqslant L_{k}\right]+\left[G_{i}<L_{k}\right]\left[O_{i} \geqslant L_{k}\right]\right)} $

(8)

式(4)—(8)中，N为样本总数，G_i为格点降水值，O_i为观测降水值，G和O分别为格点和观测降水的平均值。k代表降水分级检验等级，L_k为第k个降水等级区间的下界，[ ]代表逻辑转数值的运算符，逻辑值为正时取1，否则取0。

2 结果与分析

采用STMAS方法，基于上述设计的融合试验方案，选取江西省2019年主汛期5—6月强降水时段气象站、水文站观测降水以及雷达降水产品，进行逐小时融合分析，得到不同的融合分析结果，分别从空间分布对比、逐小时序列、站点独立性检验以及不同高度地形等4个方面对融合分析结果进行评估分析。

2.1 空间分布对比

将融合输入的各种降水数据源和融合分析结果进行可视化处理，用于分析各种融合结果的空间分布特征。为便于直观比较，选取同一降水时次数据进行对比分析。图 4给出2019年6月8日21时各观测降水空间分布。从图 4a、c中可见，雷达估测降水具有与站点观测降水一致的空间分布特征，其强降水落区、降水分布形态较为一致，但存在一定偏差。水文站观测降水与气象站观测降水(图 4a、b)相比，无论在降水形态还是降水量级上均较一致，说明气象站和水文站观测降水具有很好的一致性，这是两者可进行融合分析的前提和基础。

另外，图 5给出2019年6月8日21时方案Ⅰ、方案Ⅱ两组融合试验结果。从中看到，不同融合试验分析结果在空间分布上基本一致，既保留了雷达降水的精细化空间结构特征，又具有观测降水准确的数值信息，但两者融合分析结果在降水量级上有一定差别，因江西省水文站站点数量多(达3 476个)，且在空间分布上与气象站有差别，所以通过融合水文站数据可明显改变融合分析结果，差别较大的区域主要发生在强降水区域(图中黑色线圈所示)，说明通过融合水文站观测降水之后对强降水的改变较为明显。

2.2 逐小时序列对比

对江西省93个国家站点每个时次降水量进行累加求平均，得到2019年5—6月全省逐小时平均观测降水和融合分析降水量，生成观测降水和融合分析降水的逐小时序列对比图见图 6。从中看出，融合分析结果具有与观测降水一致的变化趋势，均能有效捕捉降水峰值和低值。融合分析降水较观测降水存在一定的低估，尤其在强降水时期低估程度更为明显。从数据融合原理看，多源资料融合分析是以当前分析格点为中心，在一定影响半径范围内搜索所有观测站点，利用这些站点的观测降水去调整背景场降水值来获取高精度的格点降水，由于降水存在很大的空间异质性，因此利用周边多个站点进行背景场降水调整通常会导致其在强降水偏弱、在弱降水偏强的现象。融合了高密度水文站降水的分析结果较仅融合气象站的分析结果的低估程度有明显减小，且在降水峰值时低估改进更明显。当平均小时降水量大于2 mm·h^-1时，融合了水文站和仅融合气象站的降水偏差分别为-0.23 mm·h^-1和-0.35 mm·h^-1，低估程度提高了77%，说明通过进一步融合水文站降水对提高研究区降水精度具有优势。

2.3 结果检验

为进一步对比分析融合试验结果，以江西省339个检验站(即93个国家站和246个区域站)逐小时观测降水数据为检验标准，对不同融合试验结果进行独立性定量检验评估。采用最邻近插值方法，将融合分析结果格点数据插值到检验站站点位置，得到站点上逐小时格点降水序列。利用站点的格点降水序列和观测降水序列，计算各融合分析结果的相关系数(COR)，均方根误差(RMSE)、平均误差(ME)、平均绝对误差(MAE)、降水评分(TS)等统计指标。同时根据降水强弱，本研究将站点小时观测降水量(R)分为5个等级，即0.1~1.9，2.0~4.9，5.0~9.9，10.0~19.9、≥20.0 mm·h^-1以上，分别对应小雨、中雨、大雨、暴雨和特大暴雨5个降水量级，再分别计算各融合分析结果在不同降水量级下的各项统计指标。此外，根据站点海拔高度，探讨降水融合分析在不同高度地形下的融合效果。

2.3.1 独立性检验结果

在站点小时尺度，对方案Ⅰ和方案Ⅱ的试验结果进行独立性检验。表 4给出两组方案的检验结果，从中看出，无论从全部样本还是从不同降水量级看，进一步融合水文站降水的分析结果的大部分统计指标均要优于仅融合气象站降水的分析结果，仅ME在小降水量级(0.1~1.9 mm·h^-1)下和总体样本要高于仅融合气象站的分析结果。两个融合试验的检验结果均表明ME在小降水量级为正偏差，而在其他降水量级为负偏差，也进一步说明融合降水在弱降水存在高估，而在强降水存在低估现象。从全部样本统计结果看，进一步融合水文站的融合试验结果，其COR较仅融合气象站的结果提高了9.2%，RMSE和MAE分别降低了14.3%和12.1%。从不同降水量级看，TS随降水量级增大而降低，进一步融合水文站的分析结果在大雨或以上量级的改进效果要优于小雨或中雨降水量级，在小雨量级两者的各项统计指标基本相当，COR分别为0.35和0.34，提高了2.9%，RMSE分别为1.22 mm·h^-1和1.29 mm·h^-1，降低了5.4%，而在大雨及以上量级，各项统计指标改进效果明显，其中RMSE在各量级下分别降低了19.1%、20.2%和23.1%，TS分别增加了17.2%、24.1%和39.5%。从独立检验结果看，进一步融合水文站观测降水可显著提高融合分析结果的准确性，尤其在强降水区域改进效果更为明显。

2.3.2 不同高度地形的检验结果

根据站点海拔高度(H，单位: m)，将上述339个检验站点划分为5个不同高度等级，即0≤H＜200、200≤H＜400、400≤H＜600、600≤H＜800以及800 m以上，以检验不同高度地形下的两种融合试验方案的降水融合效果。通过统计每个高度层级的站点数，分别为180、103、33、12、11个。表 5给出不同融合分析结果在不同站点高度的独立检验结果。

从表 5中可见，在高海拔地形下(高度在800 m以上)的降水融合效果较差，两组融合试验的各项统计指标均表现为最低，其COR仅为0.69和0.75，而RMSE高达1.74 mm·h^-1和1.60 mm·h^-1。对比同一高度两组融合试验结果表明，进一步融合水文站的分析结果在各高度层级下的各项统计指标均要优于仅融合气象站的分析结果，COR在各高度分别提高8.0%、10.6%、10.1%、6.1%和8.8%，而RMSE分别降低12.7%、16.8%、14.0%、14.6%、8.1%和14.1%，充分说明水文站降水不仅对低海拔地区的降水融合有明显改进，对高海拔地区也有显著的改进效果。江西省水文站点密度大，且较均匀分布在重点中小河流重要河段、山洪沟、存在病险隐患水库和大中型水闸、鄱阳湖区重点圩垸等区域，这也是通过进一步融合水文站降水可以在复杂地形下取得改进效果的主要原因，对气象站点在江西省复杂地形条件下的降水监测是一个有利补充。

3 结论

本文采用多重网格变分分析技术在1 km×1 km逐小时的时空尺度上开展了气象站、水文站与雷达估测降水的融合试验，利用2019年主汛期5—6月强降水资料，通过设计不同的融合试验方案，验证水文站降水对融合分析结果的影响，并从空间分布特征、时序变化、定量统计分析等方面对试验结果进行评估分析，主要结论如下：

(1) 气象站、水文站、雷达QPE降水等各种来源的观测降水数据具有相似的空间分布特征，其强降水落区、降水分布形态高度一致。不同融合试验的分析结果在空间分布上基本一致，融合了水文站与仅融合气象站降水有一定差别，差别较大的区域主要分布在强降水区域。

(2) 通过与站点降水对比，不同融合试验的分析结果具有一致的变化趋势，但均存在一定程度的低估，而进一步融合了水文站降水的低估程度明显减少。

(4) 独立检验结果表明，无论从总体样本还是在不同降水量级上，融合了水文站的分析结果的各项统计指标均要优于仅融合气象站的分析结果，相关系数提高了9.2%，均方根误差和平均绝对误差分别降低了14.3%和12.1%。

(5) 从不同高度地形的检验结果看，进一步融合了水文站降水不仅对低海拔地区的降水融合有明显改进，对高海拔地区也有显著的改进效果。相关系数在不同地形条件下分别提高了8.0%、10.6%、10.1%、6.1%和8.8%，而均方根误差分别降低了12.7%、16.8%、14.0%、14.6%、8.1%和14.1%。

综上所述，通过进一步融合水文站降水资料，可以显著提高研究区的降水融合分析的准确性，尤其对强降水的把握更加准确。本研究在已开展的气象、水文降水资料对比分析基础上，通过气象站、水文站降水资料的融合应用，发挥了高密度水文雨量站监测资料优势，对提高降水过程特别是灾害性暴雨的监测能力具有重要意义。本研究旨在探讨进一步融合水文站资料对降水融合分析的改进效果，不足之处是未对本研究结果与国内同类降水产品进行对比分析，有待于下一步收集整理相关数据产品对试验结果进行全面检验评估。

致谢：感谢国家气象信息中心多源数据融合团队在数据融合方法实现过程中给予的帮助和指导。