基于贝叶斯滤波的通用数据同化方法技术

一种基于贝叶斯滤波的通用数据同化方法：在预报步骤中将初始值集合输入分析模型获得预报集合值；在更新步骤中使用集合卡尔曼滤波计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每个预报集合；或采用粒子滤波利用集合预报值计算每个集合样本的重要性权重，利用归一化重要性权重计算有效粒子数，根据权重对集合进行重采样得到更新后的分析值和分析集合；或采用无迹卡尔曼滤波计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每个预报集合；再重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。使地球遥感观测信息和陆面过程模型信息有效融合，形成误差较小的陆面过程信息预报系统。

【技术实现步骤摘要】

本专利技术涉及地球系统科学信息处理领域，具体地说是涉及一种基于贝叶斯滤波的通用数据同化算法，使得地球遥感观测信息和陆面过程模型信息能够有效融合，从而能够形成误差较小的陆面过程信息(如土壤水分、土壤温度等)的预报系统。
技术介绍
陆面数据同化的核心思想是在陆面过程模型的动力框架内，通过数据同化的算法融合不同来源和不同分辨率的直接与间接观测，将陆面过程模型和各种观测算子(如辐射 传输模型)集成为不断地依靠观测而自动调整模型轨迹，并且减小误差的预报系统。现代数据同化算法一般可分为连续同化和顺序同化两大类。其中，连续同化是指在一个同化窗口内，利用优化算法，通过迭代而不断调整模型初始场，最终将模型轨迹拟合到离散的观测点(多维的观测向量)上。连续同化的经典算法为变分法，常用的变分同化方法有3DVar和4DVar ;顺序同化则一般利用滤波算法，在有观测的时刻，利用观测真值在误差加权的基础上对模型状态进行更新，从而获得模型状态的后验优化估计。顺序同化一般采用卡尔曼滤波和扩展卡尔曼滤波。近年来，以集合卡尔曼滤波为代表的非线性滤波方法也得到了广泛的应用，尽管国内外学者对数据同化算法进行了大量的研究工作，但目前还存在以下不足I)无论是集合卡尔曼滤波还是经典的扩展卡尔曼滤波，都假设误差的先验概率分布为多维高斯分布，因此，这两种算法在解决非高斯问题的时候精度不够好；另外，需要着重指出的是服从高斯分布的变量在经过非线性系统的转换以后也会呈现非高斯分布的特征。2)变分同化方法(例如3DVar和4DVar) —般需要求模型的伴随算子，但在陆面数据同化领域，由于大部分陆面过程模型和辐射传输模型都是非线性的，而且陆面过程模型广泛使用阈值而导致其不连续性(邱崇践，1997)，因此陆面过程模型的伴随算子在通常情况下很难得到(Reichle et al.，2002)。3)在滤波同化方法中，扩展卡尔曼滤波作为非线性滤波领域最有名的算法(Anderson & Moore, 1979)在数据同化领域得到了十分广泛的应用。但是，当这种算法应用于非线性特征明显的数据同化系统时，却遇到了很大的困难。首先EKF滤波需要对模型算子和观测算子做线性化变换，这种线性化对复杂的模型系统是非常困难的，而且，线性化只是原模型在局部的一阶近似，模型动力学的细节(高阶信息)都丢失了。另外在高维数据同化系统中对于EKF滤波，误差协方差矩阵的计算量和存储量都非常大，阻碍了 EKF滤波的应用。信号处理领域在发展非线性非高斯滤波算法方面起步早，知识积累非常丰富，有许多方法都可被陆面数据同化所借鉴。这些方法包括基于确定性采样的方法和基于随机采样的方法，其中，确定性采样的方法包括无迹卡尔曼滤波和中心差分卡尔曼滤波；随机采样的方法包括基于蒙特卡罗随机采样的采样重要性重抽样粒子滤波和无迹粒子滤波。这些新方法与集合卡尔曼滤波相比，都有一定的优势。但是，仍需要采用更精确的算法，使数据同化系统的精度不断提闻。
技术实现思路
本专利技术的目的在于克服上述已有技术的不足，从贝叶斯(Bayes)估计的角度对现有非线性滤波方法进行抽象化，从理论和主要计算方法上对他们进行统一，建立基于贝叶斯滤波的通用数据同化算法。本专利技术是这样实现的利用C++程序设计语言编程实现集合卡尔曼滤波、无迹卡尔曼滤波、和粒子滤波 等贝叶斯滤波方法，由于贝叶斯滤波涉及大量矩阵运算，程序中的矩阵运算调用高性能的数学函数库BLAS (基本线性代数子例程)和LAPACK (线性代数程序包)来实现，并且利用OpenMP发展并行算法以提高系统的效率。编程实现以马其赛特旋转(Mersenne Twister)算法为核心的长周期的高性能随机数生成器，包括单变量正态分布、多变量正态分布、偏正态分布等各种概率分布函数的简单随机数生成算法；将地统计学的随机模拟方法与简单随机数生成算法相结合来构造集合卡尔曼滤波和粒子滤波算法所需的初始样本(图I)。利用函数指针的形式统一各个算法的接口，可以为以后的数据同化研究提供一个通用、快速、稳健、能够自动纠错、物理约束更强的同化系统算法平台。本专利技术提供的具有三种算法，都是基于贝叶斯理论推导出来的。其中第一种，主要包括预报和更新两个步骤，其中预报步骤将初始值集合输入分析模型，获得预报集合值；更新步骤集合卡尔曼滤波利用预报集合计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每一个预报集合，得到更新后的分析值和分析集合；重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。所述的，其中，预报步骤具体为x! = /(x! ，v!)/tl/t-1 J v /t-1 5 /t-1Jr n A(A)rP = v ) N _ IA = I NX = — > X叫-I N叫_1更新步骤依次为权利要求1.一种，主要包括预报和更新两个步骤，其中 预报步骤将初始值集合输入分析模型，获得预报集合值； 更新步骤集合卡尔曼滤波利用预报集合计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每一个预报集合，得到更新后的分析值和分析集合； 重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。2.根据权利要求I所述的，其中，预报步骤具体为3.一种，主要包括预报和更新两个步骤，其中 预报步骤将初始值集合输入分析模型，获得预报集合值； 更新步骤粒子滤波利用集合预报值计算每个集合样本的重要性权重，利用归一化重要性权重计算有效粒子数，根据权重对集合进行重采样，减少权重小的样本，利用权重较大的样本进行补充，得到更新后的分析值和分析集合；粒子的权重的取值范围为0-1，权重越接近于0说明粒子的权重越小； 重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。4.根据权利要求3所述的，其中，预报步骤依次为 (1)初始化从初始概率密度函数P(Xtl)中采样得到样本式 (2)预报从转移概率密度函数尸( 14—0中采样得到样本4 更新步骤依次为 (3)计算重要性权重=(o'k_lP(y k\x'k) N (4)归一化重要性权重=WKM=——l-—— (5)计算有效粒子数eff)2 / 二 I (6)重采样如果小于预先给定的重采样阈值，则根据粒子的权重公丨对其重采样，得到N个新的粒子，这些粒子的概率密度函数为p (xk I yk) (7)重新设置每个粒子的权重Fori = I, . . . , N,使=(b'k =HN (8)计算滤波后的状态后验估计值Sk二^ / 二 I5.一种，主要包括预报和更新两个步骤，其中 预报步骤将初始值集合输入分析模型，获得预报集合值； 更新步骤无迹卡尔曼滤波利用预报集合计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每一个预报集合，得到更新后的分析值和分析集合； 重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。6.根据权利要求5所述的，其中，预报步骤依次为 (1)选择参数k，a和P的值。参数k(k> 0)确保协方差矩阵的半正定，一般情况下k = 0 ;参数a (0 < a < I)控制点的本文档来自技高网...

【技术保护点】
一种基于贝叶斯滤波的通用数据同化方法，主要包括预报和更新两个步骤，其中：预报步骤：将初始值集合输入分析模型，获得预报集合值；更新步骤：集合卡尔曼滤波利用预报集合计算预报误差协方差矩阵，通过观测值和卡尔曼增益矩阵更新每一个预报集合，得到更新后的分析值和分析集合；重新将更新后的分析集合作为分析模型的初始值进行下一步预测和同化，重复上述预报步骤和更新步骤。

【技术特征摘要】

【专利技术属性】
技术研发人员：韩旭军，李新，摆玉龙，
申请(专利权)人：中国科学院寒区旱区环境与工程研究所，
类型：发明
国别省市：