CN105069295A - 基于卡尔曼滤波的卫星以及地面降水测量值同化方法 - Google Patents

Abstract

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

Description

基于卡尔曼滤波的卫星以及地面降水测量值同化方法

技术领域

本发明涉及基于卡尔曼滤波的卫星以及地面降水测量值同化方法。

背景技术

在大气科学、海洋科学、环境科学的研究中，降雨量是个非常重要的物理量。降水体系在整个大气的能量传输中起着十分重要的作用，它在相当程度上引起了全球湿度和热量的重新分布，因此降水量的精确测量对进行气象灾害预报、减灾防灾和安排工农业生产等各项工作都具有重要意义，而且将推动大气海洋耦合的研究，帮助我们深入理解潜热在大尺度大气环流中的作用，进一步提高对暴雨系统的认识和分析能力。特别对于目前的全球环境气候变化问题，全球的降水分布及垂直结构的遥感更是研究中的关键环节，也是薄弱环节。目前只有少数国家可以获得可靠的全球降雨分布数据。

除了这些科学的意义外，降雨量还可以给农业、水资源管理、水灾预报带来许多实际的效益，因此精确地测量全球降雨量是安排工农业生产、进行气象灾害预报、减灾防灾等各项工作的实际需要。在与人们生活密切相关的日常天气预报中，如何提高预报的准确率是一个最受关注的问题，而降水预报则是最重要也是最困难的问题。

测量降水的方法很多，一种是利用地面雨量计。地面雨量计作为一种直接的降水测量方法，其优点是对观测点来说，具有较高的准确性。缺点是由于降水在时间和空间上存在很大的不均匀性，而雨量计分布密度有限，不能反映较大面积的降水分布，尤其对于偏远地区，盆地，以及地形复杂的区域。另一种是利用遥感卫星进行探测。遥感卫星的优点是具有较大的空间覆盖性，能通过主动和被动探测反演降水方法使其可以对全球降水分布做出时间和空间上都较为连续的观测。早前的红外/可见光主要是通过云探测对降水进行间接反演，其基本原理是依据云顶温度越低，降雨量则越大来估计雨强。但其缺点是可见光和红外波对云的穿透性较差，所获得的主要信息来自降水云顶部，降低了遥感信息与地面观测资料的可比性。简单总结，来自雨量计的数据在测量空间上是精确的，但可能在时间和数量上不精确；而雷达和卫星测量在时间上通常是精确的，但在空间上是欠精确的，在数量上更需要进行进一步的校验。

近年来，随着星载微波探测器的发展，微波资料反演降水成为卫星遥感降水量的一个重要方向，因为相比热红外通道和可见光通道而言，微波能深入一定的云层、甚至到达地表，可直接反映降水云的微物理特性，与降水的关系更为直接。涌现了很多星载卫星，包括PERSIANN(Sorooshianetal.,2000),CMORPH(Joyceetal.,2004),和PERSIANN-CCS(Hongetal.,2004)。热带降雨测量计划测雨雷达(简称TRMM)是第一个主动遥感的星载降水雷达，其多卫星降雨分析(TMPA)系统是迄今为止对降雨数据有着最优预测的系统，尤其是降雨卫星产品——TMPA准实时3B42RT。该数据集是一套覆盖南北纬60度之间的全球3小时降雨平均格点降雨数据集，水平分辨率为0.25°×0.25°时间分辨率为3小时。该算法提供了一种融合了TRMM实时合并了的被动微波方法和微波校准红外算法。目前的方法是简单的替换：对于每一个格点的值，如果HQ值有效，就用它替换；否则就用VAR值替换。

3B42RT算法的输入就是我们上面讲的融合TRMM实时3B40RT值和3B41RT值，它兼顾了两者，该方法同样适用于其他网格化的“高质量”降水估计或者其它来源的“低质量”降水估计,只要这两套输入很好地相互校准。

3B42RT是实时估计TRMM和其他数据系统的3级产品。该系统发展成为一个开发融合近似全球降水估计和利用其接近实时的有效的输入的优势的平台。

尽管TRMM-TMPA系统帮助我们获得更精准的降雨测量信息，但还是不可避免的存在系统偏差。

发明内容

针对现有技术的不足，本发明提供一种基于卡尔曼滤波的卫星以及地面降水测量值同化方法，通过卡尔曼滤波的方法对地面雨量计和TRMM星载降水雷达的实时降雨测量数据进行数据同化，从而更精准的预估出降雨量。

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

一种基于卡尔曼滤波的卫星以及地面降水测量值同化方法，包括以下步骤：

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

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

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

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

Z_t＝G_t/S_t(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)所述通过卡尔曼滤波方法对测量偏差预估值和噪声误差方差预估值进行同化具体为：

将测量偏差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_t)P_t^{-}\\ K_t=P_t^{-}{(P_t^{-}+R)}^{-1}\end{array}\right.---\left(5\right)$

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

本发明的有益效果是：

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

2、采用卡尔曼滤波和标准的TRMM3B42RT产品相比，是更有效的数据同化工具，在日分布和月分布的降雨数据中都和地面观测值有着更密切的一致性；

3、卡尔曼滤波的同化方法不仅适用于雨量站分布密度大的区域，在雨量计分布密度小的复杂地形区域也有着较强的适用性；

4、基于卡尔曼滤波的数据同化方法在秋季降雨测量中有更好的优化效果。

附图说明

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

图1是本发明基于卡尔曼滤波的数据同化方法的处理过程的流程图；

图2a是在整个泾河流域卫星测量数据与地面站点(日尺度)数据的比较图；

图2b是在整个泾河流域同化结果与地面站点(日尺度)数据的比较图；

图2c是在整个泾河流域卫星测量与地面站点(月尺度)数据的比较图；

图2d是在整个泾河流域同化结果与地面站点(月尺度)数据的比较图；

图3a是在整个泾河流域同化过程中随卡尔曼滤波参数改变的CC变化比较图；

图3b是在整个泾河流域同化过程中随卡尔曼滤波参数改变的RMSE变化比较图；

图3c是在整个泾河流域同化过程中随卡尔曼滤波参数改变的ME变化比较图；

图3d是在整个泾河流域同化过程中随卡尔曼滤波参数改变的BIAS变化比较；

图4a是划定泾河流域雨量站空间分布，在泾河流域上中下流域各取3个网格获得卫星测量数据与地面站点(日尺度)数据的比较图；

图4b是划定泾河流域雨量站空间分布，在泾河流域上中下流域各取3个网格获得同化结果与地面站点(日尺度)数据的比较图；

图4c是划定泾河流域雨量站空间分布，在泾河流域上中下流域各取3个网格获得卫星测量与地面站点(月尺度)数据的比较图；

图4d是划定泾河流域雨量站空间分布，在泾河流域上中下流域各取3个网格获得同化结果与地面站点(月尺度)数据的比较图；

图5a是划定泾河流域雨量站空间分布，泾河流域上中下流域各取3个网格获得的同化过程中随卡尔曼滤波参数改变的CC变化比较图；

图5b是划定泾河流域雨量站空间分布，泾河流域上中下流域各取3个网格获得的同化过程中随卡尔曼滤波参数改变的RMSE变化比较图；

图5c是划定泾河流域雨量站空间分布，泾河流域上中下流域各取3个网格获得的同化过程中随卡尔曼滤波参数改变的ME变化比较图；

图5d是划定泾河流域雨量站空间分布，泾河流域上中下流域各取3个网格获得的同化过程中随卡尔曼滤波参数改变的BIAS变化比较；

图6a是减少同化所用雨量站数量后，卫星测量数据与地面站点(日尺度)数据的比较图；

图6b是减少同化所用雨量站数量后，同化结果与地面站点(日尺度)数据的比较图；

图6c是减少同化所用雨量站数量后，卫星测量与地面站点(月尺度)数据的比较；

图6d是减少同化所用雨量站数量后，同化结果与地面站点(月尺度)数据的比较图；

图7a是减少同化所用雨量站数量后，随卡尔曼滤波参数改变的CC变化比较图；

图7b是是减少同化所用雨量站数量后，随卡尔曼滤波参数改变的RMSE变化比较图；

图7c是减少同化所用雨量站数量后，随卡尔曼滤波参数改变的ME变化比较图；

图7d是减少同化所用雨量站数量后，随卡尔曼滤波参数改变的BIAS变化比较图；

图8a是泾河流域从2006年1月到2008年12月卫星测量和高密度站点同化分析后的月平均降水量值；

图8b是泾河流域从2006年1月到2008年12月卫星测量和稀疏站点同化分析后的月平均降水量值。

具体实施方式

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

如图1至图8(b)所示，基于卡尔曼滤波的卫星以及地面降水测量值同化方法，包括以下步骤：

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

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

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

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

(1)预处理：

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

星载降雨雷达(TRMM3B42RT，以下简称TRMM)为数据同化提供实时的遥感测量数据。并用S_t表示。数据同化的另一来源为地面雨量计，雨量计的数据来源于中国气象机构管辖的国家标准站。因此可确保两个数据来源相互独立。网格化之后的地面雨量计测量的数据用G_t表示。

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

Z_t＝G_t/S_t(1)

(2)时间更新

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

由于现在无法确定不同时刻的TRMM的日测量偏差，假定预估值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)$

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

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

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

(3)测量更新

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

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

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

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

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

在泾河流域进行试验，泾河流域是典型的半干旱气候，位于南纬34°46’-37°19’,东经106°14’-108°42’。每年平均温度8℃，平均降雨量539.1mm,，平均地表径流18.32mm。该流域的地表雨量计分布相互独立且共有200个地面雨量计每天记录降雨数据。

为了对比分析卫星降雨测量数据的空间误差和经过同化后的结果，特选定几个典型的指标：CC(相关系数),RMSE(均方根误差),ME(平均误差),和BIAS(偏差)。两组数据的日分布特征散点图如图2(a)和图2(b)所示。由图像可得知，卫星雷达的RMSE值，ME值,和BIAS值均高于同化数据的值，但卡尔曼滤波的CC值显著提高。这种特征在月分布的特征散点图中更为显著，如图2(c)图2(d)所示。为进一步对卫星测量数据和同化数据进行对比，特在该流域的上流域，中部流域和下流域各取三个网格，获得的同化前后数据结果如图4(a)-图4(b)所示。

对比结果表明：在所选择的网格上，同化数据在日监测或月监测中都比卫星降雨测量数据校准率更高，偏差值更小，。

2、在地面雨量计分布少的地区的同化结果分析

在缺乏雨量计的地区，常常使用例如TRMM之类的气象卫星来提供大范围的降雨预测，然而，这样的方法在定点上预测不准确，且缺乏持续性。为模拟该过程，我们将使用的泾河流域的雨量计降雨测量值从200个降到13个，使用同样的同化技术，所获得的两组数据的对比图如图6(a)-图6(d)所示。

结果表明：减少地面雨量计数量对同化结果影响较小，说明基于卡尔曼滤波的同化方法具有较好的稳定性和适应性，能在雨量计分布密度小的区域发挥作用。

3、同化结果的参数敏感度分析

在进一步分析同化研究中参数的变化，选取了TRMM的测量偏差预估值噪声方差预估值(P_t),测量噪声方差(R),状态噪声方差(Q)及卡尔曼滤波增益(K_t)作为卡尔曼滤波过程中敏感参数。为了简化分析，取定中A为单位矩阵，且在仿真实现的过程中，设定这五个参数在0.1-0.9之间变化，各情况下泾河流域降雨测量数据的参数敏感性测试如图3(a)-图3(d)，图5(a)-图5(d)，图7所示，图3是在整个泾河流域，同化过程中各指标随卡尔曼滤波参数改变的变化曲线；图5是划定泾河流域雨量站空间分布，泾河流域上中下流域各取3个网格获得的同化过程中各指标随卡尔曼滤波参数改变的变化曲线；图8是减少同化所用雨量站数量后，同化过程中各指标随卡尔曼滤波参数改变的变化曲线。

通过分析可发现，在卡尔曼滤波同化的过程中，最敏感的参数是测量偏差预估值和卡尔曼滤波增益K_t，是指TRMM的测量偏差，因此在卡尔曼滤波同化过程中越大，则说明TRMM测量值和雨量计测量值之间的误差越小。K_t主要用来表示在同化过程中平衡前一时刻的预估偏差值和测量偏差值。取两种极端情况，当K_t＝1时表明预估测偏差值与测量偏差值相等则前一时刻的预估信息在同化过程中没有起作用，当K_t＝0时表明没有输入测量偏差值，因此，和K_t的取值对于同化结果有着至关重要的影响，在研究流域取定在0.4-0.7之间，在0.1-0.3之间和使得同化结果的相关系数更优，均方根误差更小。

4、不同季节的同化结果分析

图8(a)表示有着200个雨量计站点的泾河流域从2006年1月到2008年12月的月平均降水量测量值，图8(b)表示的是泾河流域雨量计站点减少为13个时的月平均降水量测量值。结果表明：(1)卫星降雨测量数据结果和地面雨量计的测量结果相比，总体偏高。(2)同化之后的数据在降雨测量数据对比中获得了最优效果，且在雨量计分布较少的区域同样适用。(3)地面雨量站的密度分布大小对同化结果影响不大，同化技术具有稳健性。

进一步的研究，对比分析雷达降雨测量结果和同化结果的季节性分布降水量(见表1)可知，同化研究方法对于秋季有着更显著的优化效果。

表1

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

Claims (4)

1.基于卡尔曼滤波的卫星以及地面降水测量值同化方法，其特征在于，包括以下步骤：

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

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

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

2.根据权利要求1所述的基于卡尔曼滤波的卫星以及地面降水测量值同化方法，其特征在于：所述步骤(1)获得相互独立的卫星测量数据和地面雨量站数据，同时对测量值比较分析得出测量偏差具体为：通过星载降雨雷达为数据同化提供实时的遥感测量数据，并将其记作S_t，将网格化之后的地面雨量计测量的数据记作G_t；将星载降雨雷达提供的实时遥感测量数据与地面雨量计的测量值之间的测量偏差记作Z_t，t表示时间，将Z_t表示为一个自回归模型如下式：

Z_t＝G_t/S_t(1)。

3.根据权利要求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是状态噪声方差。

4.根据权利要求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_t)P_t^{-}\\ K_t=P_t^{-}{(P_t^{-}+R)}^{-1}\end{array}\right.---\left(5\right)$

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