CN105975791A - 一种基于双平滑法的稀疏区域降水估算融合方法 - Google Patents

Abstract

本发明公开了一种基于双平滑法的稀疏区域降水估算融合方法，基于数据同化思想，将雨量计观测降水点数据与TMPA栅格降水数据融合在一起形成一组离散点，最终将两种融合的数据集转变为一种改进后的新数据集；同时，本发明采用平滑处理来纠正TMPA融合前的不连续性。双平滑方法在降水值和空间分布上更加合理有效，并且在雨量计分布稀疏地区平均偏差和均方根误差随着时间变化较小，估算强降水更为准确，同时降水场可视化结果也更加平滑。本发明方法尤其适用于西部地区自动气象站空间分布稀疏的地区。

Description

一种基于双平滑法的稀疏区域降水估算融合方法

技术领域

本发明涉及一种基于双平滑法的稀疏区域降水估算融合方法，属于统计数据融合方法。

背景技术

在天气变化要素中，降水是一种非常重要的气象要素。降水估算则是研究水文气象中一个极为重要的参考量。及时准确地估算降水对于国民经济、社会生活、人民生命财产安全都有重要的意义。传统降水估算是利用地面观测资料获取，观测站观测是最直接、有效的获取降水资料的方法。很多水文模式的初值都是用流域上分布的观测网的观测资料，用加权平均法，平滑函数插值法或者地理统计方法获得。事实上，我国地面观测站主要分布在东南和中部地区，其他区域的空间分布比较稀疏。由于降水空间变异性很大，目前的观测网络通过插值得到的面雨量估算有很大的不确定性。

近年来，科学家们为了能够采用覆盖到全球的卫星图像进行降水估算，进行了大量的研究。热带降水测量卫星(TRMM)，提供了丰富的全球热带雨量资料。目前，TMPA(实时3小时的TRMM多卫星降水分析)已经直接应用到许多水文气象研究中，但是，TMPA的空间分辨率对于需要获取气象或者水文异质性的研究是不理想的。卫星估算的低分辨率成为水文气象研究的限制，但因其有延时、时间采样频率和全球覆盖率高的特点，使TMPA是值得采用的。

20世纪末，空间信息统计学方法被引入到融合卫星与地面观测估算降水算法中，沈艳等人对融合产品质量评估结果表明，融合产品有效利用地面观测和卫星反演降水各自的优势，在降水量值和空间分布上均更加合理，减少了平均偏差和均方根误差，使得产品质量进一步提高，而且在定量监测强降水具有优势。

发明内容

发明目的：针对自动气象站分布稀疏地区现有降水数据获取方法存在的不足，本发明提供一种基于双平滑法的稀疏区域降水估算融合方法，是一种融合TMPA和地面观测数据的新统计方法，即使存在TMPA数据偏差和没有较强空间假设的情况下也可以表现出很好的非参数统计框架，使得在自动气象站稀疏区域能够很好地间断性校正和空间插值。相比沈艳等人的方法，本发明方法的平均偏差、均方根误差和相关系数都有明显改善，估算降水量更准确并且具有更好的视觉表现。

技术方案：为实现上述目的，本发明采用的技术方案为：

一种基于双平滑法的稀疏区域降水估算融合方法，包括如下步骤：

(1)利用TMPA(实时3小时的TRMM多卫星降水分析)的移动平均线在卫星影像上产生一个平滑场，该平滑场内包含的TMPA栅格总数为n，TMPA栅格的边长为L，第i个TMPA栅格记为S_i；移动窗口分辨率与TMPA栅格分辨率完全相同，移动窗口的边长为2L，将移动窗口划分为四个方形区域，分别记为P₁、P₂、P₃、P₄，第j个方向区域记为P_j；

移动窗口沿移动平均线移动，S_i与P_j的相交区域记为A_ij；以T(S_i)表示S_i处的TRMM3B43降水估计，则S_i处的平滑TMPA降水估算T^*(S_i)为：

$T^{*}\left(S_i\right)=\underset{j=1}{\overset{4}{\mathrm{\Σ}}}{w}_{ij}T\left(S_i\right)---\left(1\right)$

$w_{ij}=A_{ij}/\underset{t=1}{\overset{n}{\mathrm{\Σ}}}{A}_{tj}---\left(2\right)$

将S_i处的平滑TMPA降水估算T^*(S_i)称为背景场X_B(S_i)；

(2)根据背景场X_B(S_i)和基于雨量计得到的实测场X_O(S_i)，计算残差场D(S_i)：

D(S_i)＝X_B(S_i)-X_O(S_i) (3)

(3)利用双平滑技术对残差场D(S_i)进行空间插值加权平均，获得背景误差场

(4)根据背景场X_B(S_i)和背景误差场获得估算场X_M(S_i)：

$X_M\left(S_i\right)=X_B\left(S_i\right)-{\widehat{\mathrm{\μ}}}_B^D\left(S_i\right)---\left(4\right)$

估算场X_M(S_i)即为融合结果。

具体的，所述步骤(3)具体包括如下步骤：

(31)将残差场D(S_i)转换成栅格化伪数据

$\hat{D}\left(S_i^{*}\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{K}_1\left(\right||S_i^{*}-S_i||/h_1)D\left(S_i\right)}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{K}_1\left(\right||S_i^{*}-S_i||/h_1)}---\left(5\right)$

其中：为S_i对应的伪数据栅格，K₁为高斯核函数，带宽h₁为经验值(这里取h₁＝0.3)；

(32)通过扩充残差场D(S_i)和栅格化伪数据组成的数据集估算得到背景误差场

${\widehat{\mathrm{\μ}}}_B^D\left(S_i\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{K}_2\left(\right||\stackrel{\‾}{S}-S_i||/h_2)D\left(S_i\right)+\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{K}_2\left(\right||\stackrel{\‾}{S}-S_k^{*}||/h_2)\hat{D}\left(S_k^{*}\right)}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{K}_2\left(\right||\stackrel{\‾}{S}-S_i||/h_2)+\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{K}_2\left(\right||\stackrel{\‾}{S}-S_k^{*}||/h_2)}---\left(6\right)$

其中：K₂为Epanechnikov核函数，带宽h₂为经验值(这里取h₂＝0.3)；从所有n个TMPA栅格中选取出K个目标TMPA栅格总数，第k个目标TMPA栅格记为S_k，为S_k对应的伪数据栅格。

有益效果：本发明提供的基于双平滑法的稀疏区域降水估算融合方法，基于数据同化思想，将雨量计观测降水点数据与TMPA栅格降水数据融合在一起形成一组离散点，最红将两种融合的数据集转变为一种改进后的新数据集；同时，本发明采用平滑处理来纠正TMPA融合前的不连续性。双平滑方法在降水值和空间分布上更加合理有效，并且在雨量计分布稀疏地区平均偏差和均方根误差随着时间变化较小，估算强降水更为准确，同时降水场可视化结果也更加平滑。本发明方法尤其适用于西部地区自动气象站空间分布稀疏的地区。

综上所述，融合TMPA和地面观测数据至少有两个优点：(1)两种雨量估算误差是独立的，融合这两种数据可以弥补估测误差，提高精度；(2)在地面观测站分布稀疏区域，TMPA构建的降水区域具有重要的参考价值。

附图说明

图1为本发明的基于双平滑法的稀疏区域降水估算融合方法流程示意图；

图2为全国自动气象站雨量计分布图；

图3为2009年7月两种方法在C1地区的效率系数；

图4为2005到2010年5-9月沈艳等人方法求得的效率系数空间分布图；

图5为2005到2010年5-9月本发明方法求得的效率系数空间分布图；

图6为2005到2010年5-9月沈艳等人方法在稀疏地区C1的偏差空间分布图；

图7为2005到2010年5-9月本发明方法在稀疏地区C1的偏差空间分布图；

图8为2005到2010年5-9月采用沈艳等人方法生成的月均降水场栅格图；

图9为2005到2010年5-9月采用本发明方法生成的月均降水场栅格图。

具体实施方式

下面结合实施例对本发明作更进一步的说明。

本试验选用我国雨量计降水数据和TRMM卫星降水数据。首先利用年卫星降水数据TRMM 3B43平滑融合后的TMPA(实时3小时的TRMM多卫星降水分析)数据和雨量计数据由公式(3)，获得残差场D(S_i)，再由公式(5)得到栅格化伪数据然后扩充残差场D(S_i)和栅格化伪数据组成的数据集由公式(6)估算得到背景误差场最后将所得的背景场X_B(S_i)和背景误差场由公式(4)获得估算场X_M(S_i)，此过程如图1所示。

全国总共有839个自动气象站，该数据文件时间段选取2005到2010年。由于我国地形复杂，自动气象站雨量计网络分布不均，总体表现为东部密集西部稀疏的分布特征，因此需要评估不同区域站网密度下的融合方法质量。根据实验区选取三个同等大小的网格区：西部C1区有12个站点，站点分布最稀疏；其次是中部C2区，有68个站点；东南C3区站点分布最为密集，有162个站点，如图2所示。将C1、C2和C3分别与TRMM 3B43降水产品进行融合，并检验3种区域站网密度下的融合效果。

表1给出了2009年7月三种中国区域站网密度下融合方法的平均偏差、均方根误差和相关系数。随着雨量计分布密度的增加，融合产品的平均偏差和均方根误差减小，空间相关性增大。本发明方法在雨量计分布稀疏的C1地区平均偏差为-0.082mm/h，均方根误差为1.960mm/h，相关系数为0.455；在雨量计分布一般的C2地区估算误差减小，平均偏差和均方根误差分别减小到-0.058mm/h和1.833mm/h，相关系数则增加到0.701；在雨量计分布密集的C3地区，平均偏差、均方根误差进一步减小，而相关系数可以达到0.722。相比沈艳等人方法，本发明方法估算降水没有较大起伏，整体估算更加稳定。在C1地区，采用本发明方法后平均偏差由-0.121mm/h改进到-0.082mm/h，相关系数也从0.309提高到0.455，说明对于我国雨量计分布稀疏地区降水估算有显著改善。

表1 2009年7月3种不同区域(图1中的C1，C2，C3)站网密度下融合方法的平均偏差、均方根误差和相关系数

效率系数常常作为水文模型的效率评价指标。本发明利用效率系数进一步试验比较两种融合方法在稀疏地区C1的模型效率，如图3所示为2009年7月两种融合方法效率系数的比较。从图中可以看到，所有点都位于1:1线上方，这也意味着Y轴采用的方法优于X轴采用的方法，因此，试验证明在稀疏地区采用本发明的融合方法更加有效。

为了更好地说明在雨量计分布稀疏地区本发明方法估算降水更为准确，将对比分析各降水量级下的两种融合方法的统计指标，来试验不同降水量级下的融合结果。根据降水的强弱可以把小时降水量分成5个等级，小于1.0mm、1.0-2.5mm、2.5-8.0mm、8.0-16.0mm和大于16.0mm。由表2可见，沈艳等人方法的平均偏差和相对偏差由正值变为负值，其说明高估了小于1.0mm/h的低值降水，低估了大于等于1.0mm/h的降水，并且随着雨强的增大，偏差也逐渐增大。而采用本发明的融合方法后的降水场虽然还是对降水估算有所偏差，但是同一量级下的平均偏差、相对偏差和均方根误差的绝对值比沈艳等人方法明显减小。其中，采用本发明的融合方法的降水量级大于8.0mm/h的改进结果最明显，尤其是大于16mm/h的强降水，相对偏差从之前的-30.967％减小到-14.461％，对于强降水估算更为准确。

表2 2005到2010年5—9月不同降水量级下融合的统计结果

下面探讨和比较空间尺度上的模型有效性。图4和图5分别表示2005到2010年5-9月沈艳等人方法和本发明方法求得的效率系数空间分布。我国有些偏远地区，雨量计分布十分有限，例如，图4和图5中C1地区的雨量计分布较为稀疏，沈艳等人的方法在该地区不能很好地估算降水，效率系数大致在0.5以下，而采用本发明方法在该地区的效率系数范围可以达到0.5至0.9之间甚至部分地区达到0.9以上，这充分表明了本发明方法在雨量计分布稀疏地区估算降水更加准确有效。

为了进一步说明在雨量计分布稀疏地区本文融合方法得到的降水估算较为精确，利用两种方法分别计算2005到2010年5-9月两种方法所得降水估算与C1地区站点观测降水的偏差(Bias)，结果如图6和图7所示。从图6和图7可以看出，两种方法在西部稀疏地区均低估了降水，这可能与该地区雨量计分布稀疏和局部对流活动有关；采用沈艳等人的方法所得偏差范围主要集中在小于-1，仅有2个站点偏差在-1至-0.1之间；而采用本发明双平滑方法处理后的偏差范围则大致在-1至-0.1，和-0.1至0之间，分别有站点数5个和6个，还有一个偏差范围在0至0.1之间。这进一步说明了在稀疏地区双平滑融合方法可以更有效地减小误差，从而能够更好地反映降水情况。

图8和图9分别表示2005到2010年5-9月采用沈艳等人方法和本发明方法生成的月均降水场栅格图。从图中可以看出，两种方法生成的降水场整体分布很相似，发生大降水事件主要集中在东南地区；由沈艳方法产生的降水面较粗糙，特别是在地面观测站稀疏地区，如C1区，而本发明方法在稀疏地区采用双平滑处理后提高了降水场的精度，在图中表现得更加平滑。

实验结果表明，相比沈艳等人方法采用双平滑方法在降水值和空间分布上更加合理有效，并且在雨量计分布稀疏地区平均偏差和均方根误差随着时间变化较小，估算强降水更为准确，并且降水场可视化也更加平滑。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (2)

1.一种基于双平滑法的稀疏区域降水估算融合方法，其特征在于：包括如下步骤：

(1)利用TMPA的移动平均线在卫星影像上产生一个平滑场，该平滑场内包含的TMPA栅格总数为n，TMPA栅格的边长为L，第i个TMPA栅格记为S_i；移动窗口分辨率与TMPA栅格分辨率完全相同，移动窗口的边长为2L，将移动窗口划分为四个方形区域，分别记为P₁、P₂、P₃、P₄，第j个方向区域记为P_j；

移动窗口沿移动平均线移动，S_i与P_j的相交区域记为A_ij；以T(S_i)表示S_i处的TRMM 3B43降水估计，则S_i处的平滑TMPA降水估算T^*(S_i)为：

$T^{*}\left(S_i\right)=\underset{j=1}{\overset{4}{\Σ}}{w}_{ij}T\left(S_i\right)---\left(1\right)$

$w_{ij}=A_{ij}/\underset{i=1}{\overset{n}{\Σ}}{A}_{tj}---\left(2\right)$

将S_i处的平滑TMPA降水估算T^*(S_i)称为背景场X_B(S_i)；

(2)根据背景场X_B(S_i)和基于雨量计得到的实测场X_O(S_i)，计算残差场D(S_i)：

D(S_i)＝X_B(S_i)-X_O(S_i) (3)(3)利用双平滑技术对残差场D(S_i)进行空间插值加权平均，获得背景误差场

(4)根据背景场X_B(S_i)和背景误差场获得估算场X_M(S_i)：

$X_M\left(S_i\right)=X_B\left(S_i\right)-{\widehat{\mathrm{\μ}}}_B^D\left(S_i\right)---\left(4\right)$

估算场X_M(S_i)即为融合结果。

2.根据权利要求1所述的基于双平滑法的稀疏区域降水估算融合方法，其特征在于：所述步骤(3)具体包括如下步骤：

(31)将残差场D(S_i)转换成栅格化伪数据

$\hat{D}\left(S_i^{*}\right)=\frac{\underset{i=1}{\overset{n}{\Σ}}{K}_1\left(\right||S_i^{*}-S_i||/h_1)D\left(S_i\right)}{\underset{i=1}{\overset{n}{\Σ}}{K}_1\left(\right||S_i^{*}-S_i||/h_1)}---\left(5\right)$

其中：为S_i对应的伪数据栅格，K₁为高斯核函数，带宽h₁为经验值；

(32)通过扩充残差场D(S_i)和栅格化伪数据组成的数据集估算得到背景误差场

${\widehat{\mathrm{\μ}}}_B^D\left(S_i\right)=\frac{\underset{i=1}{\overset{n}{\Σ}}{K}_2\left(\right||\stackrel{\‾}{S}-S_i||/h_2)D\left(S_i\right)+\underset{k=1}{\overset{K}{\Σ}}{K}_2\left(\right||\stackrel{\‾}{S}-S_k^{*}||/h_2)\hat{D}\left(S_k^{*}\right)}{\underset{i=1}{\overset{n}{\Σ}}{K}_2\left(\right||\stackrel{\‾}{S}-S_i||/h_2)+\underset{k=1}{\overset{K}{\Σ}}{K}_2\left(\right||\stackrel{\‾}{S}-S_k^{*}||/h_2)}---\left(6\right)$

其中：K₂为Epanechnikov核函数，带宽h₂为经验值；从所有n个TMPA栅格中选取出K个目标TMPA栅格总数，第k个目标TMPA栅格记为S_k，为S_k对应的伪数据栅格。