CN105425319A - 基于地面测量数据校正的降雨卫星暴雨同化方法 - Google Patents

Abstract

本发明是基于地面测量数据校正的降雨卫星暴雨同化方法，包括以下步骤：预处理：获得相互独立的卫星测量数据和地面雨量站数据，同时对测量值比较分析得出测量偏差；时间更新：将测量值估值和噪声误差估值都表示为自相关函数；测量更新：通过卡尔曼滤波方法对测量偏差预估值和噪声误差方差预估值进行同化；本发明不仅订正了系统的误差，消除部分噪音影响，并且能在强降雨气候进行更加精准度预测，对降雨量估计值的可靠性进行了进一步的提升。

Description

基于地面测量数据校正的降雨卫星暴雨同化方法

技术领域

本发明属于降雨量预估技术领域，尤其及一种基于地面测量数据校正的降雨卫星暴雨同化方法。

背景技术

降水是地球水循环系统的重要组成部分，降水量的改变在建立水循环模型，维持生态环境的稳定，农业生产，洪水报警预测方面起着十分重要的作用。精准地测量降水及其区域和全球分布，长期以来一直是一个颇具挑战性的科学研究目标。

传统的雨量站和气象雷达系统是一种直接获取降雨数据的方式，但是由于降雨在时间和空间上的不均匀性以及测量站点的分布局限性，很难获得大区域和全球性的降雨分布，尤其对于偏远地区，盆地，以及地形复杂的区域。因此，近年来，具有全球覆盖性和几乎不间断的时间连续性的卫星降雨预估产品在弥补地面雨量计和地基雷达降雨测量的局限性方面的作用越来越受到重视。

目前主流的卫星降水产品包括红外主导的降水预估，该产品对使用人工神经网络获得的遥感信息进行降水反演，还有后来的PERSIANN-CSS(Hongetal.,2004)，微波主导的CMORPH(Joyceetal.,2004)和基于TRMM的降雨卫星产品3B42RT,3B42V6(Huffmanetal.,2007)等，这些降雨卫星产品已经在灾害监测，水文建模和气象监测方面获得了广泛应用。热带降雨测量计划测雨雷达(简称TRMMPR)是第一个主动遥感的星载降雨雷达，其多卫星降雨分析(TMPA)系统是迄今为止对降雨数据有着最优预测的系统。TMPAZ准实时3B42RT将微波辐射计放置在多个低轨道卫星上，全球降雨观测卫星能以每三小时为周期提供观测数据，其数据集是一套覆盖南北纬60度之间的全球3小时降雨平均格点降雨数据集。然而，当要求获得单位小时的降雨数据时则会产生明显的抽样误差。因此，利用填缝技术将卫星数据与来自地球同步卫星的红外数据相结合，以获得具有高时间分辨率的全球降雨图变得十分必要。

全球降水卫星制图(GSMaP；Okamotoetal.2005；Kubotaetal.2007；Aonashietal.2009；Ushioetal.2009)是近年来解决基于卫星的降雨预估在时间上和空间上具有高尺度变异性的有效方法。此技术起始于2002年，由日本科学技术振兴机构(JST)和日本宇宙航空机构(JAXA)的降水测量科学团队提供支持。此项目的主要目的是使用当前在轨的被动微波辐射计数据生产具有高精度，高分辨率的全球降雨产品；此外，改善降水物理模型和降雨速率反演算法；并进一步评估产品精度，为全球降雨测量(GMP)项目的开展奠定基础。这个项目广泛收集了来自被动微波(PMW)和红外(IR)传感器的卫星预测输入数据，可以通过将被动微波传感器的数据和红外辐射计的数据相结合的技术，提供具有较高时间精度和空间精度的产品。基于此技术的全球降雨预测产品在洪水预警系统等水文领域有广泛应用。

未经过雨量计修正的全球降雨卫星制图(GSMaP)产品叫做GSMaP_MVK，它能提供表面降雨率为0.1度，时间分辨率为一小时的全球降雨预测值。为了提高卫星降水可靠性，一种有效的方法是将此测量数据与另一数据来源相结合，如全球雨量计数据集。由于全球降雨卫星制图产品(GSMaP_MVK)预估和地面测量在时间和空间上抽样点不同，两组测量值会出现不匹配的情况，为了填补卫星预测和地面雨量站之间由时空分辨率不同引起的预测误差，基于雨量计的全球降雨卫星制图校准产品GSMaP_Gauge(Kubotaetal.，2007)应运而生。

对该产品可能产生的误差进行量化并进一步进行记录，可增大其应用价值和应用范围。然而，基于该产品现有的降水评估主要局限在日本领域(Ushioetal.2009；Kubotaetal.2009)，其他地域仍缺乏详细的研究。此外，在具有强降雨的案例中，该产品的技术方法在预测应用中比起观测站较大程度上降低了降雨预估值。

发明内容

针对现有技术的不足，为了水文学的进一步研究及工程应用，数据同化是修正调整实时测量数据的一个有效方法。因此，本文的主要目的就是给出一个方法，就是将全球降雨卫星制图产品的预测值和一定密度的雨量计观测所预测的降雨量数值进行同化，从而更精准的预估出降雨量，尤其是极端强降雨的情况。

本发明解决其技术问题所采用的技术方案如下，

一种基于地面测量数据校正的降雨卫星暴雨同化方法，包括以下步骤：

(1)预处理：获得相互独立的卫星测量数据和地面雨量站数据，同时对测量值比较分析得出测量偏差，具体为：GSMaP_MVK为数据同化提供实时的遥感测量数据，并将其记作S_t，将网格化之后的地面雨量计测量的数据记作G_t，将GSMaP_MVK提供的实时遥感测量数据与地面雨量计的测量值之间的测量偏差记作Z_t，将Z_t表示为一个自回归模型如下式：

(2)时间更新：将测量值估值和噪声误差估值都表示为自相关函数；

(3)测量更新：通过卡尔曼滤波方法对测量偏差预估值和噪声误差方差预估值进行同化。

进一步的，所述步骤(2)中将测量值估值和噪声误差估值都表示为自相关函数具体为：

设Z_t的预估值为x_t，并且x_t满足自相关函数关系：

x_t＝Ax_t-1+B_t(2)

其中，A是测量偏差的一阶相关系数，B_t是一个均值为0的独立随机变量；将x_t的预估值表达为：

${\hat{x}}_t=A{\hat{x}}_{t-1}---\left(3\right)$

将B_t预估方差P_t ^-表示为：

P_t ^-＝AP_t-1A^T+Q(4)

其中，P_t-1是t-1时刻预测的误差方差，Q是状态噪声方差。

进一步的，步骤(3)所述通过卡尔曼滤波方法对测量偏差预估值和噪声误差方差预估值进行同化具体为：

将测量偏差估计值用表示，噪声误差估计值用P_t表示，使用卡尔曼滤波更新方法得到以下表达式：

其中，K_t表示t时刻的卡尔曼滤波增益，R表示测量噪声方差。

本发明的有益效果是：

1、能在获取卫星实时数据的基础上结合雨量站信息获得更可靠的降雨数据，同时减小数据误差；

2、卡尔曼滤波后的同化数据和标准的全球降水卫星制图产品GSMaP_MVK相比，是更有效的数据同化工具，在降雨强度较大的情况中和地面观测值有着更密切的一致性。

3、卡尔曼滤波的同化方法不仅适用于海拔比较低的地区，在海拔较高的山地也同样有着较强的适用性；

4、基于卡尔曼滤波的数据同化方法在捕捉连续时间序列的降水时空变化中有更好的优化效果。

附图说明：

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明基于地面测量数据校正的降雨卫星暴雨同化方法的处理过程的结构示意图；

图2(a)是整个泾河流域在0.1°*0.1°的分辨率下，a,b,c类网格区域内卫星测量值随累计降水量和海拔高度变化的CC参数变化三维图；

图2(b)是整个泾河流域在0.1°*0.1°的分辨率下，a,b,c类网格区域内同化结果随累计降水量和海拔高度变化的CC参数变化三维图；

图2(c)是整个泾河流域在0.1°*0.1°的分辨率下，a,b,c类网格区域内卫星测量值随累计降水量和海拔高度变化的RMSE参数变化三维图；

图2(d)是整个泾河流域在0.1°*0.1°的分辨率下，a,b,c类网格区域内同化结果随累计降水量和海拔高度变化的RMSE参数变化三维图；

图2(e)是整个泾河流域在0.1°*0.1°的分辨率下，a,b,c类网格区域内卫星测量值随累计降水量和海拔高度变化的BIAS参数变化三维图；

图2(f)是整个泾河流域在0.1°*0.1°的分辨率下，a,b,c类网格区域内同化结果随累计降水量和海拔高度变化的BIAS参数变化三维图；

图3(a)是在整个泾河流域从2010年07月22日早上6点到2010年07月24日中午11点，a类网格区域(即累计降水量超过100mm的区域)的雨量计降雨数据，卫星降雨数据和同化后的降雨数据之间的对比图；

图3(b)是在整个泾河流域从2010年07月22日早上6点到2010年07月24日中午11点，b类网格区域(即累计降水量小于100mm且同化后的数据偏差有改善的区域)的雨量计降雨数据，卫星降雨数据和同化后的降雨数据之间的对比图；

图3(c)是在整个泾河流域从2010年07月22日早上6点到2010年07月24日中午11点，c类网格区域(即累计降水量小于100mm且同化后的数据偏差更大的区域)的雨量计降雨数据，卫星降雨数据和同化后的降雨数据之间的对比图；

具体实施方式

下面结合附图对本发明的优选实施例进行详细阐述，以使本发明的优点和特征能更易于被本领域技术人员理解，从而对本发明的保护范围做出更为清楚明确的界定

如图1至图3(c)所示，基于地面测量数据校正的降雨卫星暴雨同化方法，包括以下步骤：

(1)预处理：获得相互独立的卫星测量数据和地面雨量站数据，同时对测量值比较分析得出测量偏差。

(2)时间更新：由于卡尔曼滤波方程是递推形式，且仅与前一时刻的有关信息相关，因此将测量值估值和噪声误差估值都表示为自相关函数。

(3)测量更新：通过卡尔曼滤波方法对测量偏差和噪声值进行同化。

下面对这三个步骤进行详细说明：

(1)预处理：

卡尔曼滤波校准法是对一个变量的两个独立估计进行加权平均。当选择合适的权重因子使得平均结果具有最小方差时，平均结果就是卡尔曼滤波器的输出。

全球降水卫星制图产品(GSMaP_MVK)为数据同化提供实时的遥感测量数据。并用S_t表示。数据同化的另一来源为地面雨量计，雨量计的数据来源于中国气象机构管辖的国家标准站，因此可确保两个数据来源相互独立。将来自网格化之后的地面雨量计测量的数据用G_t表示。

GSMaP_MVK的降雨预估值S_t与相应的雨量计预估值G_t之间的测量偏差用Z_t表示。如果我们假定GSMaP_MVK降雨预测偏差在空间上保持不变，且遵循同一个马尔科夫进程，则测量偏差Z_t将是一个自回归模型，其参数将由卡尔曼滤波器不断进行更新。

(2)时间更新

卡尔曼滤波方程的解是递推形式，而且仅与前一时刻的有关信息相关。因此，想要获得t时刻的偏差滤波输出值，需要首先设定当前时刻的测量偏差值Z_t和对应点误差方差。

由于现在无法确定GSMaP_MVK的测量偏差的时间变化规律，先假定预估值x_t满足自相关函数关系

x_t＝Ax_t-1+B_t(2)

其中A是测量偏差的一阶相关系数，B_t是一个均值为0的独立随机变量，即高斯白噪声。因此x的预估可表达为：

${\hat{x}}_t=A{\hat{x}}_{t-1}---\left(3\right)$

此时(3)式提供了计算预估值的基本方法，预估值不仅考虑了新的观察测量值，也能纠正在新的测量值中可能出现的任何错误。

B_t预估方差可参照误差方差，用P_t ^-表示：

P_t ^-＝AP_t-1A^T+Q(4)

其中，P_t-1是t-1时刻预测的误差方差，Q是状态噪声方差。

(3)测量更新

测量更新即根据当前时刻t的实际测量值，用卡尔曼滤波的方法对测量偏差预估值和噪声误差方差的预估值P_t ^-进行同化研究。更新完成后的测量偏差估计值用表示，噪声误差估计值用P_t表示。由于x的真实值和测量值之间存在测量误差，用z_t表示t时刻检测到的测量偏差。使用卡尔曼滤波的测量更新之后将得到以下表达式：

其中，K_t表示t时刻的卡尔曼滤波增益，R表示测量误差方差，表示(z_t-x_t)的是时变方差。

下面通过几个方面来对本方法的可靠性进行验证。

1、在泾河流域的同化结果分析

在泾河流域进行试验，泾河流域是典型的半干旱气候，位于南纬34°46’-37°19’,东经106°14’-108°42’。每年平均温度8℃，平均降雨量539.1mm,，平均地表径流18.32mm。该流域海拔范围从海平面350米的通道出口到超过2900米上游山区，流域地形从西北到东南显著下降。该流域的地表雨量计分布相互独立且共有200个地面雨量计每天记录降雨数据。

在研究过程中，为了更好地进行对比研究，根据雨量站所在的位置将划分的网格分为三类，a类是累计降雨量超过100mm的网格点；b类是累计降水量小于100mm，但同化后的偏差统计数据有改善的网格点；c类是累积降水量小于100mm,且同化后的统计数据偏差更大的网格点。

将GSMaP_MVK，GSMaP_Gauge,同化处理之后的降雨数据结果和地面雨量计的测量结果进行对比(如表1(a)-表1(c))。

总体而言你，对于全球降水卫星制图产品GSMaP_MVK，位于东部区域的网格点的测量数据准确度优于西部区域。同类网格点(如a类)之间的对比表明，在任何区域的网格点，全球降水卫星制图产品GSMaP_Gauge和同化过程处理过的数据都有了极大的改善，大大降低了偏差率和RMSE(均方根误差)。

对比结果表明，用雨量计的数据对卫星降水数据进行同化处理在降雨量较大的案例中能有效提高全球降水卫星制图的精确性。另一方便，通过比较b类网格和c类网格可得知，在降水量较少时，全球降水卫星制图产品GSMaP_MVK所获得的数据的精确度较高，即GSMaP_MVK数据的精确率将随着降水量的增加而降低。

由表3可知，采用同化技术的产品和全球降水卫星制图产品GSMaP_MVK相比都低估了降雨量，采用同化技术获取的数据偏差为-14.03％，GSMaP_MVK的偏差为-56.03％，两者的相关系数也相差无几，但是采用同化技术获得的数据的均方根误差参数RMSE较之得以显著改善；与具有内置校准的全球降水卫星制图产品GSMaP_Gauge相比，两者在CC,RMSE,BIAS等参数指标中性能相近，都有了大大改善。

2、对不同海拔区域的同化结果分析

为了进一步进行性能比较，特选定几个典型的指标：CC(相关系数),RMSE(均方根误差),和BIAS(偏差)，如表表2(a1)-表2(c3)所示。

在海拔较高的地区的降水更值得注意，因为海拔越高其降水量通常越大。表2(a1)-表2(c3)中，a类网格点的平均海拔为1729米，平均降水量为175.41mm，b类网格点的平均海拔为1336.33米，平均降水量为75.36mm，c类网格点的平均海拔为1271.67米，平均降水量为47.42mm。此外，高海拔地区的降水量往往在预测时偏低。表2(a1)-表2(c3)中全球降水卫星制图产品GSMaP_MVK的平均降水偏差在三类网格中分别为-77.53％,-27％和-6.01％。在同类网格点(如a类)的对比中，同样能发现这一特点，海拔较高的15和12号网格的统计结果比海拔较低的17号网格的偏差率高。

结果表明，全球降水卫星制图产品GSMaP_MVK在海拔较低的区域能获得较为精准的数据，但随着海拔的升高，数据精准率将降低。

3、对不同降水量的情形下的同化结果分析

为进一步研究，按照累计降水量和海拔的不同分别对三个降水预估指标(CC、RMSE、BIAS)分析。通过对比可以发现，卫星降雨预测在高海拔的强降雨地区不能取得较好的结果，而在海拔较低，降水量较少的地区能较好的预测。同时在图2中可以发现，在a类网格点中，同化后的预测数据大大减小了偏差值。

4、在连续降水时间序列中的同化结果分析

图3描述了泾河流域2010年7月22日早上6点到2010年7月24日中午11点的连续降水情况。图3(a)是a类网格点中，雨量计，卫星和同化方法三者之间各参数的对比图，图3(b)是b类网格点中，雨量计，卫星和同化方法三者之间各参数的对比图，图3(c)是c类网格点中，雨量计，卫星和同化方法三者之间各参数的对比图。

对比结果表明，(1)和地面雨量站的测量数据相比，卫星产品GSMaP_MVK的降雨预测结果整体偏低，将降雨量增大时更为突出。(2)就相关系数而言，同化处理后的数据具有更高的精准度，包括在预测条件不利的c类网格点内，如图3(c)。(3)通过空间气象站同化技术得出的降雨数据在暴雨的情况下能显著提升卫星降雨预测水平，在a类网格点内即平均降水量超过100mm的区域，表现更为显著。

表1(a)

表1(b)

表1(c)

表2(a1)

表2(a2)

表2(a3)

表2(b1)

表2(b2)

表2(b3)

表2(c1)

表2(c2)

表2(c3)

表3

表1(a)是划定泾河流域雨量站空间分布，a类网格点区域内不同海拔的雨量计降雨数据，卫星降雨数据，带校准的卫星降雨数据和同化结果之间的累计降水量数值比较；

表1(b)是划定泾河流域雨量站空间分布，b类网格点区域内不同海拔的雨量计降雨数据，卫星降雨数据，带校准的卫星降雨数据和同化结果之间的累计降水量数值比较；

表1(c)是划定泾河流域雨量站空间分布，c类网格点区域内不同海拔的雨量计降雨数据，卫星降雨数据，带校准的卫星降雨数据和同化结果之间的累计降水量数值比较；

表2(a1)是划定泾河流域雨量站空间分布，a类网格点内卫星降雨数据，带校准的卫星降雨数据和同化结果随站点海拔变化的CC值变化的对比；

表2(a2)是划定泾河流域雨量站空间分布，a类网格点内卫星降雨数据，带校准的卫星降雨数据和同化结果随站点海拔变化的RMSE值变化的对比；

表2(a3)是划定泾河流域雨量站空间分布，a类网格点内卫星降雨数据，带校准的卫星降雨数据和同化结果随站点海拔变化的BIAS值变化的对比；

表2(b1)是划定泾河流域雨量站空间分布，b类网格点内卫星降雨数据，带校准的卫星降雨数据和同化结果随站点海拔变化的CC值变化的对比；

表2(b2)是划定泾河流域雨量站空间分布，b类网格点内卫星降雨数据，带校准的卫星降雨数据和同化结果随站点海拔变化的RMSE值变化的对比；

表2(b3)是划定泾河流域雨量站空间分布，b类网格点内卫星降雨数据，带校准的卫星降雨数据和同化结果随站点海拔变化的BIAS值变化的对比；

表2(c1)是划定泾河流域雨量站空间分布，c类网格点内卫星降雨数据，带校准的卫星降雨数据和同化结果随站点海拔变化的CC值变化的对比；

表2(c2)是划定泾河流域雨量站空间分布，c类网格点内卫星降雨数据，带校准的卫星降雨数据和同化结果随站点海拔变化的RMSE值变化的对比；

表2(c3)是划定泾河流域雨量站空间分布，c类网格点内卫星降雨数据，带校准的卫星降雨数据和同化结果随站点海拔变化的BIAS值变化的对比；

表3是在整个泾河流域，卫星降雨数据，带校准的卫星降雨数据和同化结果的平均降水预估指标的对比。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何不经过创造性劳动想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书所限定的保护范围为准。

Claims (3)

1.基于地面测量数据校正的降雨卫星暴雨同化方法，其特征在于，包括以下步骤：

(1)预处理：获得相互独立的卫星测量数据和地面雨量站数据，同时对测量值比较分析得出测量偏差，具体为：

通过全球降水卫星制图产品GSMaP_MVK为数据同化提供实时的遥感测量数据，并将其记作S_t，将网格化之后的地面雨量计测量的数据记作G_t；将GSMaP_MVK提供的实时遥感测量数据与地面雨量计的测量值之间的测量偏差记作Z_t，t表示时间，将Z_t表示为一个自回归模型如下式：

$Z_t={\log }_{10}\left(\frac{G_t}{S_t}\right)---\left(1\right);$

(2)时间更新：将测量值估值和噪声误差估值都表示为自相关函数；

(3)测量更新：通过卡尔曼滤波方法对测量偏差预估值和噪声误差方差预估值进行同化。

2.根据权利要求1所述的基于地面测量数据校正的降雨卫星暴雨同化方法，其特征在于，所述步骤(2)中将测量值估值和噪声误差估值都表示为自相关函数具体为：

设Z_t的预估值为x_t，并且x_t满足自相关函数关系：

x_t＝Ax_t-1+B_t(2)

其中，A是测量偏差的一阶相关系数，B_t是一个均值为0的独立随机变量；

将x_t的预估值表达为：

${\hat{x}}_t=A{\hat{x}}_{t-1}---\left(3\right)$

将B_t预估方差表示为：

$P_t^{-}={\mathrm{AP}}_{t-1}{A}^T+Q---\left(4\right)$

其中，P_t-1是t-1时刻预测的误差方差，Q是状态噪声方差。

3.根据权利要求1所述的基于地面测量数据校正的降雨卫星暴雨同化方法，其特征在于，步骤(3)所述通过卡尔曼滤波方法对测量偏差预估值和噪声误差方差预估值进行同化具体为：

将测量偏差Z_t的预估值用表示，噪声误差预估值用P_t表示，使用卡尔曼滤波更新方法得到以下表达式:

$\left\{\begin{array}{c}{\hat{x}}_t={\hat{x}}_t^{-}+K_t(z_t-{\hat{x}}_t^{-})\\ P_t=(1-K)P_t^{-}\\ K_t=P_t^{-}{(P_t^{-}+R)}^{-1}\end{array}\right.---\left(5\right)$

其中，K_t表示t时刻的卡尔曼滤波增益，R表示测量噪声方差。