一种结合最优观测时间窗口的强耦合数据同化方法与流程

本发明属于耦合气候模式系统的强耦合数据同化与数值预报技术领域，具体涉及一种结合最优观测时间窗口的强耦合数据同化方法。

背景技术：

耦合气候模式能合理模拟大气，海洋，陆地以及海冰等地球系统中各个子系统之间的相互作用过程，从而模拟气候变化过程。但耦合模式中又往往存在各种误差与不确定性(例如不完美的数值实现和物理参数化过程以及不合理的模式内部参数值等)，从而造成模式输出结果偏离真实的气候特征和变化过程的实际观测数据，从而限制耦合模式的气候预报能力。传统弱耦合数据同化过程在各个子系统独立进行状态估计(子系统内的观测信息只用来估计自身的状态)基础上利用模式的耦合动力学特征在子系统之间传输观测信息，从而保证模式各个子系统能得到一致和连贯的调整。传统强耦合数据同化过程则在弱耦合同化过程的基础上，利用海洋慢变属性(当前时刻的海洋状态中包含之前时刻的大气状态信息，所以当前时刻海洋状态与之前时刻的大气状态之间相关性更大)这一特点，取之前时刻的大气观测均值来实现对海洋状态的进一步的估计与优化，维持耦合模式内部时间尺度交互的特性。缺陷在于受计算条件所限，同化周期远远大于观测时间间隔，弱耦合同化过程中只考虑同化时间点上的观测信息，忽略了非同化时间点上的观测信息，造成了观测信息的极大浪费，限制了弱耦合同化过程中的耦合模式状态估计的精度；从而影响了对弱耦合同化过程中大气与海洋状态相关性的估计与采样精度，进而限制了强耦合同化过程中大气观测对海洋状态的估计精度，限制了强耦合同化过程的同化质量以及后续耦合模式的大气与海洋的数值预报能力。

设立观测时间窗口通常用于收集有效的观测信息以用于每一个数据同化与周期，从而提高每一个模式状态估计周期内的观测利用率。一方面传统方法中将观测时间窗口内所有观测均权重化能在减小计算代价的同时减小观测的误差的影响，缺陷是观测时间窗口内观测与状态估计时刻的模式背景场存在时间差异性，且观测对应模式状态与背景场之间的相关性系数存在差异，传统的均值化会引入上述时间差异与相关性差异，从而造成观测时间窗口内观测信息被低估或者高估，从而降低耦合模式状态的准确性。另一方面观测时间窗口的大小与同化时间间隔大小之间的关系直接决定了背景场与观测之间的相关性，传统方法中扰动观测又会引入观测扰动误差，如何权衡两者之间的影响也是引入观测时间窗口时面临的又一难题。

技术实现要素：

本发明的目的在于提供一种结合最优观测时间窗口的强耦合数据同化方法，针对传统耦合数据同化与参数优化方法中存在的观测资源利用以及状态估计与参数优化精度的不足，结合耦合模式自由积分状态的谱分析获取耦合模式状态的特征时间尺度，并依据耦合模式状态的时间尺度设立最优观测时间窗口的大小。在基于最优观测时间窗口获取有效的大气与海洋观测数据的基础上引入观测窗口内时间权重系数，最大程度上提取有效的观测信息以拟合耦合模式状态的特征变率；同时忽略模式内部参数的时变特征并引入时间窗口内的时间平均系数，实现对模式参数的更加精确的估计与优化；从而强化耦合模式的大气与海洋的数值预报能力。

本发明的目的是这样实现的：

一种结合最优观测时间窗口的强耦合数据同化方法，包括如下步骤：

步骤1：获取耦合模式状态的特征时间尺度；在耦合模式自由积分的基础上获取耦合模式状态的时间序列，并对状态的时间序列进行谱分解，以获取不同耦合状态的特征时间尺度，选取3％的特征时间尺度为对应的耦合状态设置最优观测时间窗口大小；

步骤2：获取耦合模式状态的自相关系数；在耦合模式自由积分的基础上获取耦合模式状态的时间序列，并对状态的时间序列进行自相关分析，以分别获取不同状态的自相关系数；

步骤3：对最优观测时间窗口内观测数据进行预处理；剔除最优观测时间窗口内的无效数据和异常数据，将有效的观测数据转化成同化过程所需的格式；

步骤4：构建基于最优观测时间的弱耦合同化方案，弱耦合同化方案中大气与海洋状态分别进行同化；若大气最优观测时间窗口小于大气状态的同化时间间隔，采用方案一集合平方根滤波方法，在不对观测数据进行扰动的情况下进行状态估计，否则采用方案三完整的集合卡尔曼滤波方法，对模式状态进行估计和更新；若海洋最优观测时间窗口小于海洋状态的同化时间间隔，采用方案二集合平方根滤波方法，消除观测集合扰动带来的扰动误差，否则采用方案四完整的集合卡尔曼滤波方法，对模式状态进行估计和更新；

步骤5：在模式状态数据到达稳态后计算大气与海洋状态的时间平均耦合相关系数，引入耦合相关性权重系数消除传统强耦合同化方案中差异性所带来的影响，采用方案五基于耦合相关性权重系数的强耦合同化方案，将大气观测扰动成观测集合并采用c-enkf实现利用大气观测数据对海洋状态的更新；

步骤6：采用rtp方案来对更新后的大气与海洋状态进行膨胀；

步骤7：膨胀之后的模式状态分析场为下一步模式积分与预报提供初始场，模式自由积分并进入下一个同化周期。

所述步骤4中方案一为采用集合平方根滤波方法，在不对观测数据进行扰动的情况下进行状态估计；最优观测时间窗口中的有效的大气观测数据为其中y^a,o表征大气观测，m为窗口内有效观测数量，窗口内大气状态预报场为当前时刻的大气状态模式预报场为x^a,f，观测对应大气状态自相关系数分别是(l1,…,lm)；根据观测对应状态自相关系数引入观测窗口内相关性权重系数wii为

平均化后对应的观测误差为

窗口内对应的大气状态模式预报场为为大气观测误差标准差，在集合平方根滤波方法中，模式背景场与分析场表征为集合平均与集合扰动(x)′的和：

根据集合平方根滤波方法分别对分析场集合平均和集合扰动进行更新

分析场的分析误差(1-we)p^f为

融入大气最优观测时间窗口中的所有有效观测并结合相关性权重系数得到：

大气状态分析场误差协方差为

将集合平均更新权重we与集合扰动更新权重带入更新方程对模式大气状态的集合平均与集合扰动进行更新；集合平均更新权重we与集合扰动更新权重分别为

所述步骤4中方案二为针对每一次同化周期没有海洋观测被复用，周期内模式预报数据与窗口内的观测数据之间没有相关性，此时模式背景场与观测场之间不相关；采用集合平方根滤波方法消除观测集合扰动带来的扰动误差；

最优时间观测窗口中的有效的海洋观测数据对应的海洋状态自相关系数分别是(l1,…,lm)，引入观测窗口内相关性权重系数wii，分别对模式海洋状态的集合平均与集合扰动进行更新；海洋状态x^o,f是慢变量，同化窗保持不变则最优海洋观测窗口内的海洋状态更新公式变为

x^o,a为海洋状态分析场，当前情形下的集合平均更新权重we与集合扰动更新权重分别为

为海洋状态分析场误差协方差，海洋状态更新方程为

带入更新方程中分别对模式海洋状态的集合平均与集合扰动进行更新。

所述步骤4中方案三为针对每一次同化周期都存在大气观测被复用，周期内模式预报数据与窗口内的大气观测数据之间存在相关性，此时模式背景场与观测场相关；对大气观测值进行扰动来计算背景场与观测场的相关性系数，观测扰动成集合的标准差为观测误差，采用完整的集合卡尔曼滤波方法来对模式状态进行估计和更新；

最优时间观测窗口中的有效的大气观测数据窗口内大气状态预报场为当前时刻的大气状态模式预报场为x^a,f，大气观测对应状态自相关系数分别是(l1,…,lm)，引入观测窗口内相关性权重系数wii，则扰动之后的大气观测集合为此时大气状态分析场为

推导卡尔曼增益

其中表征首先将观测扰动成观测集合，之后引入相关性权重系数平均化成一个观测数据集合，再计算平均化观测集合与背景场之间的相关性，并对模式大气状态进行更新。

所述步骤4中方案四为针对每一次同化周期都存在海洋观测被复用，周期内模式预报数据与窗口内的海洋观测数据之间存在相关性，采用完整的集合卡尔曼滤波方法对模式状态进行估计和更新；最优时间观测窗口中的有效大气观测数据扰动之后的观测集合为海洋状态分析场

推导卡尔曼增益

为海洋状态分析场误差协方差，根据上式对模式海洋状态进行更新。

所述步骤5中方案五为假设当前时刻的海洋预报场为x^o,f，窗口内大气状态预报场为强耦合同化窗口内有效的大气观测数据为观测对应状态自相关系数分别是引入观测窗口内相关性权重系数

平均化后对应的观测误差为

扰动后的大气观测集合为此时海洋状态分析场

若强耦合同化时间窗口<强耦合同化时间间隔且大气最优观测时间窗口<大气状态的同化时间间隔时，不存在大气观测被重复同化，大气预报场与观测集合之间都不存在相关性，卡尔曼增益为

否则大气预报场与观测集合之间都存在相关性，卡尔曼增益为

所述步骤6中rtp方案为：

其中与表征膨胀前模式分析场与背景场，代表膨胀后的模式分析场，a为膨胀系数。

所述步骤7中进入下一个同化周期时首先判断是否进入大气状态数据同化周期，然后判断是否进入强耦合同化周期，最后判断是否进入海洋同化周期，依此顺序依次对大气与海洋状态进行估计。

本发明有益效果在于：

(1)本发明针对传统耦合数据同化方法中存在的观测资源利用以及强耦合同化方法中平均化引入的时间差异性以及耦合相关性的差异性等不足，结合耦合模式自由积分状态的谱分析获取耦合模式状态的特征时间尺度，并依据耦合模式状态的时间尺度设立最优观测时间窗口的大小；

(2)本发明在基于最优观测时间窗口获取有效的大气与海洋观测数据的基础上，根据大气与海洋状态之间的自相关系数引入观测窗口内相关性权重系数，最大程度上提取有效的观测信息的同时减小时间与相关性差异；

(3)本发明在基于最优观测时间窗口的弱耦合同化方案的稳态大气与海洋的耦合相关性系数，在传统强耦合同化方案中引入耦合相关性权重系数，实现对耦合模式状态的更加精确的估计；从而强化耦合模式的大气与海洋的数值预报能力。

附图说明

图1为大气与海洋自相关性系数简图；

图2为大气与海洋耦合相关性系数简图；

图3为一种结合最优观测时间窗口的强耦合数据同化方法的执行流程图。

具体实施方式

下面结合附图对本发明做更进一步描述。

本发明提出一种结合最优观测时间窗口的强耦合数据同化方法，在最大程度提取观测信息的同时消除时间差异性和相关性差异性所带来的影响，以提高观测利用率的同时提高耦合模式状态估计精度，从而提高耦合模式的大气与海洋的数值预报精度。针对传统耦合数据同化方法中存在的观测资源利用以及强耦合同化方法中平均化引入的时间差异性以及耦合相关性的差异性等不足，提出一种结合最优观测时间窗口的强耦合数据同化方法。该方法首先结合耦合模式自由积分状态的谱分析获取耦合模式状态的特征时间尺度，并依据耦合模式状态的时间尺度设立最优观测时间窗口的大小；并在基于最优观测时间窗口获取有效的大气与海洋观测数据的基础上，根据大气与海洋状态之间的自相关系数引入观测窗口内相关性权重系数，最大程度上提取有效的观测信息的同时减小时间与相关性差异；同时在基于最优观测时间窗口的弱耦合同化方案的稳态大气与海洋的耦合相关性系数，在传统强耦合同化方案中引入耦合相关性权重系数，实现对耦合模式状态的更加精确的估计；从而强化耦合模式的大气与海洋的数值预报能力。

相较于传统的耦合数据同化过程，本发明最大的区别在于依据模式状态的特征时间尺度引入相应大小的最优观测时间窗口来最大程度利用观测信息；同时根据大气与海洋状态之间的自相关系数引入观测窗口内相关性权重系数以减小时间与相关性差异；同时在基于最优观测时间窗口的弱耦合同化方案的稳态大气与海洋的耦合相关性系数来更准确地估计同化过程中大气与海洋状态之间的耦合相关性，引入耦合相关性权重系数来消除传统强耦合同化方案中均值化所带来的时间与相关性差异的影响，以提高耦合模式状态估计精度。本发明一方面可以大大提高观测的利用效率，从而最大程度融入观测信息；另一方面减小窗口内观测之间的时间与相关性差异所带来的影响，实现更加精确的耦合模式状态估计；优化的耦合模式状态能为耦合模式的大气与海洋的数值预报提供更加精确的初始条件。

附图1：大气与海洋自相关系数简图(附图中黑点表征观测时间窗口内有效观测对应的自相关系数)；附图2：大气与海洋耦合相关性系数简图(附图中黑点表征强耦合窗口内有效大气观测对应的耦合相关系数，横轴向左表征大气观测时间点在当前海洋状态时间点之前)；附图3：一种结合最优观测时间窗口的强耦合数据同化方法的执行流程图。

图1与图2分别是大气与海洋自相关系数与耦合相关系数简图。下面将按照流程详述本发明提出的技术方案的具体实施方式，执行流程如图3所示。该实施方式主要包含以下几个关键内容：

步骤1)获取耦合模式状态的特征时间尺度：在耦合模式自由积分的基础上获取耦合模式状态的时间序列，并对状态的时间序列进行谱分解，以获取不同耦合状态的特征时间尺度，选取3％的特征时间尺度为对应的耦合状态设置最优观测时间窗口大小。

步骤2)获取耦合模式状态(大气与海洋)的自相关系数：在耦合模式自由积分的基础上获取耦合模式状态的时间序列，并对状态的时间序列进行自相关分析，以分别获取不同状态(大气与海洋)的自相关系数，如附图1。

步骤3)对最优观测时间窗口内观测数据进行预处理：为便于实测数据的同化，首先需要剔除最优观测时间窗口内的无效数据和异常数据，并将有效的观测数据转化成同化过程所需的格式。

步骤4)构建基于最优观测时间的弱耦合同化方案：弱耦合同化方案中大气与海洋状态分别进行同化。根据大气与海洋状态的最优观测时间窗口大小与同化时间间隔大小的对比关系，可以将基于最优观测时间的弱耦合同化方案的划分成以下四个方案：

方案1:当大气最优观测时间窗口<大气状态的同化时间间隔

当大气最优观测时间窗口<大气状态的同化时间间隔时，此时针对每一次同化周期没有大气观测被复用，周期内模式预报数据与窗口内的观测数据之间没有相关性，所以此时模式背景场与观测场之间不相关；考虑到传统enkf方法中对观测数据扰动生成观测集合过程中存在观测扰动误差，这里采用集合平方根滤波方法(ensqf)，在不对观测数据进行扰动的情况下进行状态估计。

这里假设最优观测时间窗口中的有效的大气观测数据(y^a,o表征大气观测，m为窗口内有效观测数量)，窗口内大气状态预报场为同时当前时刻的大气状态模式预报场为x^a,f，观测对应大气状态自相关系数分别是(l1,…,lm)。根据观测对应状态自相关系数引入观测窗口内相关性权重系数wii，一方面消除观测时间点与同化时间点的时间差异性，另一方面引入平均化能有效减小观测误差。则最优观测时间窗口内相关性权重系数wii为：

则平均化后对应的观测误差为：

为大气观测误差标准差。窗口内对应的大气状态模式预报场为在ensqf中，模式背景场与分析场表征为集合平均与集合扰动(x)′的和，即：

根据ensqf分别对分析场集合平均和集合扰动进行更新：

则分析场的分析误差(1-we)p^f为

所以集合平均更新权重we与集合扰动更新权重为：

这里，我们融入大气最优观测时间窗口中的所有有效观测并结合相关性权重系数，得到：

推导当前情形下的集合平均更新权重we:

大气状态分析场误差协方差为

因为观测数据与模式背景场之间没有相关性，所以集合平均更新权重we与集合扰动更新权重分别为：

带入更新方程中分别对模式大气状态的集合平均与集合扰动进行更新。

方案2:当海洋最优观测时间窗口<海洋状态的同化时间间隔

当海洋最优观测时间窗口<海洋状态的同化时间间隔时，与大气观测窗口类似，此时针对每一次同化周期没有海洋观测被复用，周期内模式预报数据与窗口内的观测数据之间没有相关性，所以此时模式背景场与观测场之间不相关；同理采用ensqf消除观测集合扰动带来的扰动误差。

假设最优时间观测窗口中的有效的海洋观测数据对应的海洋状态自相关系数分别是(l1,…,lm)，引入观测窗口内相关性权重系数wii，分别对模式海洋状态的集合平均与集合扰动进行更新。区别在于海洋状态x^o,f是慢变量，所以可以近似认为同化窗保持不变，即则最优海洋观测窗口内的海洋状态更新公式变为(x^o,a为海洋状态分析场)：

推导当前情形下的集合平均更新权重we与集合扰动更新权重分别为(为海洋状态分析场误差协方差):

则海洋状态更新方程为：

带入更新方程中分别对模式海洋状态的集合平均与集合扰动进行更新。

方案3:当大气最优观测时间窗口≥大气状态的同化时间间隔

当大气最优观测时间窗口≥大气状态的同化时间间隔时，此时针对每一次同化周期都存在大气观测被复用，周期内模式预报数据与窗口内的大气观测数据之间存在相关性，所以此时模式背景场与观测场相关，需要对大气观测值进行扰动来计算背景场与观测场的相关性系数，观测扰动成集合的标准差为观测误差。此时集合平方根滤波方法(ensqf)不再适用，这里我们采用完整的集合卡尔曼滤波方法(c-enkf)来对模式状态进行估计和更新。尽管对观测数据进行扰动会引入扰动误差，但相较于背景场与观测场之间的相关性来说影响较小。

同理假设最优时间观测窗口中的有效的大气观测数据窗口内大气状态预报场为同时当前时刻的大气状态模式预报场为x^a,f，大气观测对应状态自相关系数分别是(l1,…,lm)，引入观测窗口内相关性权重系数wii，则扰动之后的大气观测集合为此时大气状态分析场为：

推导卡尔曼增益(为大气状态分析场误差协方差)：

其中表征首先将观测扰动成观测集合，之后引入相关性权重系数平均化成一个观测数据集合。再计算平均化观测集合与背景场之间的相关性，并对模式大气状态进行更新。

方案4:当海洋最优观测时间窗口≥海洋状态的同化时间间隔

当海洋最优观测时间窗口≥海洋状态的同化时间间隔时，同理针对每一次同化周期都存在海洋观测被复用，周期内模式预报数据与窗口内的海洋观测数据之间存在相关性，并采用完整的集合卡尔曼滤波方法(c-enkf)来对模式状态进行估计和更新。同样区别在于海洋状态x^o,f是慢变量，所以可以近似认为同化窗口内海洋状态预报场保持不变，即假设最优时间观测窗口中的有效大气观测数据则扰动之后的观测集合为则海洋状态分析场

推导卡尔曼增益(为海洋状态分析场误差协方差)：

再根据上式对模式海洋状态进行更新。

步骤5)基于耦合相关性权重系数的强耦合同化方案(方案5):

在使用上述基于最优观测时间窗口的大气与海洋数据同化方案分别对海洋与大气状态进行更新，并在模式状态数据到达稳态后计算大气与海洋状态的时间平均耦合相关系数，如附图2所示。传统的强耦合同化方案中，由于海洋是慢变量以及模式大气与海洋的耦合动力学特征，当前时刻海洋状态信息中包含之前时刻的大气状态的信息(即之前时刻的大气状态与当前时刻的海洋状态的耦合相关性大于当前时刻大气与海洋状态之间的耦合相关性)。设立强耦合同化时间窗口，同时取窗口内之前时刻大气观测的均值来对当前时刻的海洋状态进行更新。之前时刻大气观测的均值较之前时刻单一大气观测具有更大的耦合相关性以及更小的观测误差。然而取之前时刻大气观测的均值忽略了观测之间的时间差异性以及耦合相关性的差异性。这里我们引入耦合相关性权重系数则能有效消除传统强耦合同化方案中上述差异性所带来的影响。

由于耦合模式中大气与海洋状态之间的耦合动力学特征，所以当之前大气状态同化了之前时刻的大气观测的同时，会将大气观测信息传递到当前时刻的海洋状态预报场中，所以无论强耦合同化频率大小，当前时刻的海洋状态预报场与强耦合同化窗口内的大气观测存在相关性，所以需要将大气观测扰动成观测集合并采用c-enkf来实现利用大气观测数据对海洋状态的更新。

假设当前时刻的海洋预报场为x^o,f，窗口内大气状态预报场为强耦合同化窗口内有效的大气观测数据为观测对应状态自相关系数分别是引入观测窗口内相关性权重系数则

则平均化后对应的观测误差为：

扰动后的大气观测集合为此时海洋状态分析场

同理推到卡尔曼增益

如果强耦合同化时间窗口<强耦合同化时间间隔且大气最优观测时间窗口<大气状态的同化时间间隔时，不存在大气观测被重复同化，则大气预报场与观测集合之间都不存在相关性，即：

否则大气预报场与观测集合之间都存在相关性，即

一般情况下强耦合同化时间间隔会介于大气与海洋状态同化时间间隔之间，而最优的强耦合同化窗口的大小也需要实验反复调整。

步骤6)为了消除强耦合同化方案中耦合协方差采样时引入的强噪声所带来的影响，我们采用rtp方案来对更新后的大气与海洋状态进行膨胀，rtp方案如下：

其中与表征膨胀前模式分析场与背景场，代表膨胀后的模式分析场，α为膨胀系数。

步骤7)膨胀之后的模式状态分析场为下一步模式积分与预报提供初始场，模式自由积分并进入下一个同化周期(首先判断是否进入大气状态数据同化周期，然后判断是否进入强耦合同化周期，最后判断是否进入海洋同化周期，依此顺序依次对大气与海洋状态进行估计)。

基于耦合状态的特征时间尺度设置相应大小的弱耦合过程最优观测时间窗口大小，在最大程度提取有效观测信息的同时实现对弱耦合过程耦合模式状态特征变率的拟合。

获取弱耦合过程耦合模式状态(大气与海洋)的自相关系数，为后续弱耦合同化过程引入观测窗口内相关性权重系数，一方面消除观测时间点与同化时间点的时间差异性，另一方面引入平均化能有效减小观测误差。

步骤4在弱耦合过程根据最优观测时间窗口与同化时间间隔之间的大小关系选取不同的同化方案，充分考虑模式背景场与观测之间的相关性以及观测扰动所带来的误差影响，实现弱耦合过程中更加精确的模式状态估计。

步骤5在弱耦合过程的基础上获取耦合状态之间的耦合相关性系数，为后续强耦合同化过程引入观测窗口内相关性权重系数，在消除观测时间点与同化时间点的时间差异性的同时有效减小观测误差；充分利用大气观测信息实现强耦合过程中对海洋状态更加精确的模式状态估计。