基于类比的离线集合卡尔曼滤波古气候数据同化系统、方法及电子设备

1.本发明涉及一种古气候场重建方法，是一种将基于类比与集合卡尔曼滤波的同化系统、方法及电子设备。


背景技术：

2.理解古气候现象能够提升对未来气候变化的预报技巧，古气候场的重建使用古气候代用指标用于评估气候状态。多数古气候场重建方法运用多元线性回归得到气候场变量，但是受限制于代用指标和气候变量间的非线性关系。
3.相应地，古气候数据同化是一种结合气候模式动力系统限制和代用观测信息，找到气候系统状态的最优估计。其中，传统的基于集合的数据同化技术使用前一个时刻通过观测更新的分析场集合作为集合模式的初始条件向前积分得到下一个循环同化的先验集合，被称为
‘
在线’循环集合同化。但是由于模式的可预报性短于观测的时间尺度，同时集合模式积分需要耗费巨大的计算量。一种
‘
离线’集合数据同化方法被提出，通过已存在的气候模拟中构造出背景场集合。
‘
离线’与
‘
在线’被系统地比较过，
‘
在线’方法并未发现比
‘
离线’方法有显著的提升。
4.因此离线的基于集合的古气候数据同化方法被广泛应用于古气候场重建中，先验集合成员常常是从已存在的长时间气候模拟中随机抓取，同时为了使重建的气候场不被每次抓取的随机性影响，相同的抓取集合在每一次同化中被重复使用。另一种基于集合的古气候数据同化方法是使用退化的粒子滤波，从集合模式的模拟中选择最优的模拟，但是即使是具备较大的集合数目，被选择的粒子也会受到维度诅咒的影响。


技术实现要素：

5.本发明针对现有技术的不足，提供了一种使用基于类比与集合卡尔曼滤波的同化方法(以下简称为aoenkf)。aoenkf先在一组长时间气候模拟中选择先验集合成员，选择标准是集合场与观测之间的均方根误差(rmse)最小或者相关系数最大。得到的先验集合成员因此捕获了
‘
流依赖’的背景场误差协方差矩阵，再将气候代用资料同化进入先验集合成员中，即可得到后验集合成员。通过本发明所述的aoenkf得到的先验集合成员，相比于
‘
在线’循环同化来说，不需要集合模式积分即可得到
‘
流依赖’的背景场误差协方差矩阵，计算量更小；相比于传统的
‘
离线’基于集合同化方法的静态的先验集合成员来说，能够更好地将代用资料信息传播给气候场。由此可知，本发明所述的同化方法能结合类比方法与集合卡尔曼滤波的优势，更好地进行古气候场的重建。
6.为实现上述技术目的，本发明将采取如下技术方案：
7.一种基于类比的离线集合卡尔曼滤波古气候数据同化方法，先根据特定标准，利用观测y在一组控制实验模拟的状态变量样本x中筛选先验集合成员，再根据集合平方根滤波的更新公式，将得到的先验集合成员通过卡尔曼增益矩阵更新后验集合平均与后验集
合扰动具体包括如下步骤：
8.步骤一、给定一组控制实验模拟的状态变量样本x及观测y
9.步骤二、在同化观测前根据不同标准选择先验集合成员
10.在同化观测前，先将状态变量样本x插值到观测y中，然后依据状态变量样本x的每一个状态变量样本xj相对于观测y的均方根误差最小或者相关系数最大的原则，从状态变量样本x中筛选出来的最初的n个按序排列的样本，构成先验集合成员；
11.步骤三、用得到的先验集合成员同化观测
12.使用平方根滤波同化误差协方差矩阵为r的观测y，先更新后验集合平均以及第i个后验集合扰动i＝1,
…
,n；然后将后验集合平均与后验集合扰动相加即可得到n个后验集合成员。
13.优选地，当从状态变量样本x中选择先验集合成员的标准为状态变量样本x 的每一个状态变量样本xj相对于观测y的均方根误差最小时，每一个状态变量样本xj相对于y的均方根误差rmse表示为：
[0014][0015]
其中：观测y表示p
×
1的观测变量；ym表示观测y中第m个观测变量；h(x
j,m
) 表示将状态变量样本xj转换到观测空间ym，即为由状态变量样本xj获得的对观测ym的估计。
[0016]
同时，n个先验集合成员的各个样本按照对应的均方根误差rmse升序排序。
[0017]
优选地，当从状态变量样本x中选择先验集合成员的标准为状态变量样本x 的每一个状态变量样本xj相对于观测y的相关系数最大时，每一个状态变量样本xj相对于y的相关系数表示为：
[0018][0019]
其中：观测y表示p
×
1的观测变量；y
′
表示观测y减去其空间平均得到的观测空间异常场；y
′m表示观测空间异常场y
′
中第m个变量；h(x
j,m
)
′
表示由状态变量样本xj获得的对观测ym估计减去其空间平均得到的观测估计值的异常场。
[0020]
同时，n个先验集合成员的各个样本按照对应的相关系数降序进行排序。
[0021]
后验集合平均的更新公式为：
[0022][0023]
后验集合扰动的更新公式为：
[0024][0025]
式中：表示步骤二所获得的n个集合先验成员的集合平均；ρ是局地化矩阵，表示舒尔积；
[0026]
xf是背景协方差pf＝xfx
ft
的平方根，的平方根，表示步骤二所得到的n个集合先验成员的第i个成员的集合扰动，h2是将模式变量映射为观测变量的观测算子，h1是观测对模式变量求偏导的雅各比矩阵。
[0027]
本发明的另一个技术目的是提供一种基于类比的离线集合卡尔曼滤波古气候数据同化系统，包括：
[0028]
获取模块，用于获取需要同化的误差协方差矩阵为r的观测y；
[0029]
同化模块，基于集合卡尔曼滤波的同化框架而搭建；
[0030]
所述同化模块在同化观测前，先将状态变量样本x插值到观测y中，然后依据状态变量样本x的每一个状态变量样本xj相对于观测y的均方根误差最小或者相关系数最大的原则，从状态变量样本x中筛选出来的最初的n个按序排列的样本，构成先验集合成员；
[0031]
所述同化模块在同化时刻，使用平方根滤波同化误差协方差矩阵为r的观测 y，先更新后验集合平均以及第i个后验集合扰动i＝1,
…
,n；然后将后验集合平均与后验集合扰动相加即可得到n个后验集合成员。
[0032]
优选地，所述同化模块从状态变量样本x中选择先验集合成员的标准为状态变量样本x的每一个状态变量样本xj相对于观测y的相关系数降序排序；每一个状态变量样本xj相对于y的相关系数表示为：
[0033][0034]
其中：观测y表示p
×
1的观测变量；y
′
表示观测y减去空间平均得到的观测空间异常场；y
′m表示观测空间异常场y
′
中第m个变量；h(x
j,m
)
′
表示由状态变量样本xj获得的对观测ym估计减去其空间平均得到的观测估计值的异常场。
[0035]
优选地，所述同化模块从状态变量样本x中选择出的n个先验集合成员的各个样本，按照每一个状态变量样本xj相对于y的均方根误差rmse升序排序；其中，每一个状态变量样本xj相对于y的均方根误差rmse为：
[0036][0037]
其中：观测y表示p
×
1的观测变量；ym表示观测y中第m个观测变量；h(x
j,m
) 表示将状态变量样本xj转换到观测空间ym，即为由状态变量样本xj获得的对观测ym的估计。
[0038]
优选地，所述同化模块在同化时刻的后验集合平均以及后验集合扰动的更
新公式分别为：
[0039][0040][0041]
式中：表示步骤二所获得的n个集合先验成员的集合平均；ρ是局地化矩阵，表示舒尔积；
[0042]
xf是背景协方差pf＝xfx
ft
的平方根，的平方根，表示步骤二所得到的n个集合先验成员的第i个成员的集合扰动，h2是将模式变量映射为观测变量的观测算子，h1是观测对模式变量求偏导的雅各比矩阵。
[0043]
本发明的第三个技术目的是提供一种电子设备，包括：至少一个处理器、至少一个存储器、通信接口和总线；其中，所述处理器、存储器、通信接口通过所述总线完成相互间的通信；所述通信接口用于该电子设备和其它电子设备的通信设备之间的信息传输；所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行上述的方法。
[0044]
根据上述的技术方案，相对于现有技术，本发明具有如下的有益效果：
[0045]
本发明基于类比的方法，通过一定的标准(相对于观测的rmse或者相关系数)选择集合先验成员，代替传统的随机抓取的静态集合先验成员来更新集合均值和集合扰动。相对于
‘
在线’循环同化的方法减少了集合模式向前积分的巨大计算量。与通常的
‘
离线’重复利用随机抓取的先验集合的同化方法相比，本发明能够构造更精确的先验集合平均，同时捕获
‘
流依赖’的背景场误差协方差信息，能够进一步减小同化后的误差。
附图说明
[0046]
图1为本发明的流程示意图。
[0047]
图2为
‘
在线’循环集合卡尔曼滤波(cenkf)、
‘
离线’集合卡尔曼滤波 (oenkf)、基于rmse标准类比的离线集合卡尔曼滤波(aoenkf_e)与基于相关系数标准类比的离线集合卡尔曼滤波(aoenkf_c)方法rmse随不同局地化值变化的序列：圆圈表示oenkf与其他试验的rmse显著不同(通过了99％的置信检验)，叉号表示aoenkf_e与aoenkf_c的rmse的显著不同(通过了99％的置信检验)。
[0048]
图3是本发明所述电子设备的结构示意图。
具体实施方式
[0049]
为了更好的理解本发明的技术方案，下面结合附图对本发明实施例进行详细描述。
[0050]
应当明确，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基
于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
[0051]
在本发明实施例中使用的术语是仅仅出于描述特定实施例的目的，而非旨在限制本发明。在本发明实施例和所附权利要求书中所使用的单数形式的“一种”、“所述”和“该”也旨在包括多数形式，除非上下文清楚地表示其他含义。
[0052]
应当理解，本文中使用的术语“和/或”仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，甲和/或乙，可以表示：单独存在甲，同时存在甲和乙，单独存在乙这三种情况。另外，本文中字符“/”，一般表示前后关联对象是一种“或”的关系。
[0053]
实施例1
[0054]
如图1所示，本发明所述的基于类比的离线集合卡尔曼滤波古气候数据同化方法，先利用观测在一组控制实验模拟的状态变量样本中根据不同的标准(模式与观测间rmse最小或者相关系数最大)筛选先验集合成员，再用得到的先验集合成员应用集合平方根滤波进行集合平均与集合扰动的更新。应用本同化方法可以进一步改进分析场。具体步骤如下：
[0055]
步骤一、给定一组控制实验模拟的状态变量及观测
[0056]
1.1.给定一组控制实验模拟的状态变量样本
[0057]
xf来自于一段控制实验模拟的状态变量中的样本xn×
t
＝{x1,
…
,x
t
},其中，t 是状态变量样本的个数；n是状态变量的维数。控制实验是由单个模式积分一段时间后得到。
[0058]
1.2.给定需要同化的观测
[0059]
给定误差协方差矩阵为r的观测y，一般假设观测误差互不相关，r是一个对角矩阵，对角线元素为观测变量的误差方差。
[0060]
步骤二、在同化观测前根据不同标准选择先验集合成员
[0061]
在同化观测前，先将状态变量插值到观测空间，然后基于最小均方根误差 (rmse)的标准或者基于状态变量与观测之间的相关系数最大的原则构建先验集合成员具体分别阐述如下：
[0062]
2.1.基于最小均方根误差标准选择先验集合成员(简称为aoenkf_e)
[0063]
aoenkf_e是选择先验集合成员基于最小均方根误差(rmse)的标准(先将状态变量插值到观测空间)。对于每一个状态变量样本xj(j＝1,
…
,t),xj相对于y的rmse可表示为
[0064][0065]
通过升序排序rmse,最初的n个样本被选为此次同化时间里的先验集合成员
[0066]
2.2.基于最大相关系数标准选择先验集合成员(简称为aoenkf_c)
[0067]
aoenkf_c选择集合先验成员基于状态变量与观测之间的相关系数(先将状态变量插值到观测空间)。对于每一个状态变量样本xj(j＝1,
…
,t),xj与y的相关系数可表示为
[0068]
[0069]
其中撇号表示通过减去空间平均得到的空间异常场。通过将相关系数进行降序排序，最初的n个样本被选为此次同化时间里的先验集合成员
[0070]
2.3.选择合适的局地化方案
[0071]
选择合适的局地化函数来局地化观测的影响，减少观测与状态变量之间的虚假相关，常用的局地化函数如gaspari和cohn(gc；gaspari和cohn 1999) 函数，由单个特征尺度参数确定。局地化函数在同化时如步骤(3.1)所示应用到卡尔曼增益矩阵上。
[0072]
步骤三、用得到的先验集合成员同化观测
[0073]
使用平方根滤波同化误差协方差矩阵为r的观测y，先更新后验集合成员的集合平均以及第i(i＝1,
…
,n)个后验集合成员的扰动，然后将后验集合平均与后验集合扰动相加得到n个后验集合成员，具体计算公式如下：
[0074]
3.1.更新集合平均
[0075]
后验集合成员的集合平均的更新公式为：
[0076][0077]
其中：ρ是局地化矩阵，表示舒尔积；表示步骤二所述的aoenkf_e或者 aoenkf_c选择得到的n个集合先验成员的集合平均；
[0078]
xf是背景协方差pf＝xfx
ft
的平方根，的平方根，表示步骤二所述的aoenkf_e或者aoenkf_c选择得到的n个集合先验成员的第i个成员的集合扰动，h2是将模式变量映射为观测变量的观测算子，h1是观测对模式变量求偏导的雅各比矩阵。
[0079]
3.2.更新集合扰动
[0080]
第i(i＝1,
…
,n)个后验集合成员的扰动更新公式为：
[0081][0082]
其中：表示步骤二所述的aoenkf_e或者aoenkf_c选择得到的n个集合先验成员的第i个成员的集合扰动；xf是背景协方差pf＝xfx
ft
的平方根，h1是观测对模式变量求偏导的雅各比矩阵；ρ是局地化矩阵，表示舒尔积。
[0083]
3.3.得到后验集合成员
[0084]
后验集合平均与后验集合扰动相加得到n个后验集合成员。
[0085]
3.4.评估同化实验结果
[0086]
当同化实验结束后得到均方根误差的时间序列进行平均作为标准评估同化实验的结果。
[0087]
实施例2
[0088]
本发明基于类比的离线集合卡尔曼滤波古气候数据同化方法同化观测，以 lorenz(2005)模型为例，以单尺度模式ii无模式误差测试了本发明的表现，并与使用传统
静态抓取，
‘
离线’集合卡尔曼滤波(oenkf，hakim 2016)方法的误差结果进行了比较。敏感试验结果表明，本发明在不同的集合大小、局地化尺度、观测误差、观测密度都优于传统的同化方法。
[0089]
步骤一、给定一组控制实验模拟的状态变量及观测
[0090]
l05模型单尺度模式ii只包含一个大尺度慢过程变量。设x为慢过程变量，单尺度模式ii可写为：
[0091][0092]
下标n代表格点序号，k为常数，f为强迫项。平流项写为：
[0093][0094]
其中∑'为一个特殊的求和算子，除首项与末项除以2外，与一般求和算子相同。 k为偶数时j＝k/2，k为奇数时j＝(k-1)/2，
[0095]
总格点数n＝960，常数k选择32。强迫项f为16。真值与控制实验的模式参数相同，表示控制实验中无模式误差存在。
[0096]
1.1.给定一组控制实验模拟的状态变量样本
[0097]
xf来自于一段控制实验模拟的状态变量中的样本,控制实验是由上述l05单尺度模式ii积分一段时间后得到。将控制实验得到的每6000积分时步(～120h) 进行平均，得到独立先验集合成员样本xn×
t
＝{x1,...,x
t
}。
[0098]
1.2.给定需要同化的观测
[0099]
给定误差协方差矩阵为r的观测y，一般假设观测误差互不相关，r是一个对角矩阵，对角线元素为观测变量的误差方差。通过向一段时间真值的平均(～120h)添加服从正态分布n(0,r)的随机扰动来制造观测。默认的观测误差方差大小为4.0。默认的观测网络为每8个格点一个观测(共120个观测格点)。每6000个积分时步的平均(～120h)产生一个观测。
[0100]
步骤二、在同化观测前根据不同标准选择先验集合成员
[0101]
2.1.基于最小均方根误差标准选择先验集合成员
[0102]
aoenkf_e是选择先验集合成员基于最小均方根误差(rmse)的标准(先将状态变量插值到观测空间)。对于每一个状态变量样本xj(j＝1,
…
,t),xj相对于y的rmse可表示为
[0103][0104]
通过升序排序rmse,最初的n个样本被选为此次同化时间里的先验集合成员
[0105]
2.2.基于最大相关系数标准选择先验集合成员
[0106]
aoenkf_c选择集合先验成员基于状态变量与观测之间的相关系数(先将状态变量插值到观测空间)。对于每一个状态变量样本xj(j＝1,
…
,t),xj与y的相关系数可表示为
[0107][0108]
其中撇号表示通过减去空间平均得到的空间异常场。通过将相关系数进行降序排序，最初的n个样本被选为此次同化时间里的先验集合成员
[0109]
2.3.选择合适的局地化方案
[0110]
gaspari和cohn(gc；gaspari和cohn 1999)函数被用作局地化函数，由单个特征尺度参数确定。局地化函数在同化时如步骤(3.1)所示应用到卡尔曼增益矩阵上。gc函数的特征尺度需要选择合适的值以获得最好的结果。
[0111]
步骤三、用得到的先验集合成员同化观测
[0112]
已有通过aoenkf_e或者aoenkf_c选择得到了n个集合先验成员以及观测y，为集合平均，是第i个成员的集合扰动，是背景协方差 pf＝xfx
ft
的平方根。观测算子h是将模式变量映射为观测变量，h是观测对模式变量求偏导的雅各比矩阵，使用平方根滤波同化误差协方差矩阵为r的观测y。
[0113]
3.1.更新集合平均
[0114]
后验集合平均的更新公式为
[0115][0116]
其中ρ是局地化矩阵，表示舒尔积
[0117]
3.2.更新集合扰动
[0118]
第i(i＝1,
…
,n)个集合成员的扰动更新公式为
[0119][0120]
3.3.得到后验集合成员
[0121]
后验集合平均与后验集合扰动相加得到n个后验集合成员。
[0122]
3.4.评估同化实验结果
[0123]
当同化实验结束后得到均方根误差的时间序列进行平均作为标准评估同化实验的结果。每种数据同化方法进行了1000月的同化实验，图2显示了各同化方法的rmse1000个月时间序列平均随不同局地化值的变化。cenkf与oenkf有着几乎相同的误差，这证明了单尺度l05模式对于初始条件几乎没有任和记忆性，模式的可预报性相比于观测的时间尺度短。而对于aoenkfs，aoenkf_e与 aoenkf_c具有类似的表现，rmse相比于oenkf显著减小。这证明了在先验集合成员中包含的
‘
流依赖’的背景场误差协方差信息能够提升同化的效果。
[0124]
实施例3
[0125]
基于上述的古气候数据同化方法，本发明提供一种基于类比的离线集合卡尔曼滤波古气候数据同化系统，包括：获取模块，用于获取需要同化的误差协方差矩阵为r的观测y；同化模块，基于集合卡尔曼滤波的同化框架而搭建；所述同化模块在同化观测前，先将状态变量样本x插值到观测y中，然后依据状态变量样本x的每一个状态变量样本xj相对于观测y的均方根误差最小或者相关系数最大的原则，从状态变量样本x中筛选出来的最初的n个按序排列的样本，构成先验集合成员；所述同化模块在同化时刻，使用平方根滤波同化误差协方差矩阵为r的观测y，先更新后验集合平均以及第i个后验集合扰动然后将后验集合平均与后验集合扰动相加即可得到n个后验集合成员。
[0126]
同化模块能够实现本发明古气候数据同化方法实施例中能够实现的各个过程，以及达到相同的有益效果，为避免重复，在此不再赘述。
[0127]
实施例四
[0128]
基于上述的古气候数据同化方法，本发明提供一种电子设备，包括：至少一个处理器、至少一个存储器、通信接口和总线；其中，所述处理器、存储器、通信接口通过所述总线完成相互间的通信；所述通信接口用于该电子设备和其它电子设备的通信设备之间的信息传输；所述存储器存储有可被所述处理器执行的程序指令，所述处理器调用所述程序指令能够执行上述的古气候数据同化方法。
[0129]
具体地，处理器用于获取需要同化的误差协方差矩阵为r的观测y；
[0130]
处理器内具有基于集合卡尔曼滤波的同化框架而搭建的同化模块；同化模块基于特定标准，利用所获取的观测y在一组控制实验模拟的状态变量样本x中筛选先验集合成员，再根据集合平方根滤波的更新公式，将得到的先验集合成员通过卡尔曼增益矩阵更新后验集合平均与后验集合扰动所述的特定标准为：依据状态变量样本x的每一个状态变量样本xj相对于观测y的均方根误差最小或者相关系数最大的原则。
[0131]
本实施例中，处理器内的同化模块，同样能够实现本发明古气候数据同化方法实施例中能够实现的各个过程，以及达到相同的有益效果，为避免重复，在此不再赘述。
[0132]
上面结合附图对本技术的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本技术的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，均属于本发明的保护之内。