CN105046046A - 一种集合卡尔曼滤波局地化方法 - Google Patents

Abstract

本发明属于实测海洋环境数据同化领域，具体涉及一种集合卡尔曼滤波局地化方法。本发明包括：对大气海洋环境实测数据进行预处理；对每个观测资料，计算观测资料处的先验观测集合成员；计算观测的先验集合平均和方差；计算集合平均的观测增量；计算各个集合扰动的观测增量；将集合平均的观测增量投影到模式状态的集合平均上；将各个集合扰动的观测增量投影到模式状态的相应集合扰动上；获得集合成员分析场。本发明对集合卡尔曼滤波中的传统局地化方法进行了改进，有效地考虑了集合平均和集合扰动所代表的不同空间尺度，显著提高了集合卡尔曼滤波的同化精度。

Description

一种集合卡尔曼滤波局地化方法

技术领域

本发明属于实测海洋环境数据同化领域，具体涉及一种集合卡尔曼滤波局地化方法。

背景技术

集合卡尔曼滤波和四维变分是目前国际上公认的两类先进的数据同化方法，它们有各自的优缺点。集合卡尔曼滤波相对于变分方法的最大优点是它通过集合抽样来模拟模式状态变量的先验概率密度分布函数，根据集合样本计算出的背景误差协方差矩阵带有模式的动力信息，因此是流依赖的。由于计算机硬件资源的限制，目前只能采用较少的集合样本(10²量级)。而对实际的海洋数值模式来说，状态变量的维数为10⁷，因此较少的集合样本就会带来显著的抽样误差。

本发明主要是对集合卡尔曼滤波中传统的局地化方案进行改进，从而提出一种新的局地化方案。根据前述，有限的集合样本会给集合卡尔曼滤波数据同化引入显著的抽样误差，比如模式状态先验方差的低估以及模式状态与远距离观测之间的伪相关性等。解决这个问题的一个有效途径是在同化过程中引入局地化机制，即一个观测仅允许影响周围一定范围内的模式格点，且距离越远影响越小。传统的集合卡尔曼滤波局地化方案将集合平均以及集合扰动看做一体进行局地化。由于集合平均和集合扰动的空间尺度存在明显的差别，传统的局地化方案存在明显的局限性。

本发明提出了一种新的集合卡尔曼滤波局地化方法，即对集合平均和集合扰动利用不同的局地化因子进行局地化，且给出了两个局地化因子的确定性关系。该方法能够大大提高集合卡尔曼滤波中传统的局地化方案的同化效果，对实现海洋环境状态场的实时构建以及预测具有重要的实用价值。

发明内容

本发明的目的在于提出一种集合卡尔曼滤波局地化方法。

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

(1)对大气海洋环境实测数据进行预处理

剔出实测数据中的无效数据和异常数据，然后将数据转化为同化所需的格式；

(2)对每个观测资料，计算观测资料处的先验观测集合成员：

对模式状态的背景集合成员，利用线性插值的方法投影到观测位置，得到观测的先验集合成员；

(3)计算观测的先验集合平均和方差：

利用步骤二得到的观测先验集合成员，计算其先验集合平均和先验方差；

(4)计算集合平均的观测增量：

利用计算出的观测的先验集合平均和先验方差、观测值以及观测误差标准差，计算集合平均的观测增量；

(5)计算各个集合扰动的观测增量：

利用计算出的观测的先验集合平均和先验方差、先验集合成员以及观测误差标准差，计算各个集合扰动的观测增量:；

(6)将集合平均的观测增量投影到模式状态的集合平均上：

利用:集合平均的局地化因子，将集合平均的观测增量线性投影到模式状态的集合平均上，完成对集合平均的同化；

(7)将各个集合扰动的观测增量投影到模式状态的相应集合扰动上：

利用集合扰动的局地化因子，将各个集合扰动的观测增量线性投影到模式状态的相应的集合扰动上，完成对集合扰动的同化；

(8)获得集合成员分析场

将由步骤(7)获得的集合扰动的分析场叠加到由步骤(6)得到集合平均的分析场上，得到最终的集合成员的分析场。

所述的步骤(4)包括：

根据计算出的观测的先验集合平均和先验标准差观测值y^o以及观测误差标准差r，按照下式计算集合平均的观测增量

$\mathrm{\Δ}\stackrel{\‾}{y}=\frac{{\left({\mathrm{\σ}}_y^p\right)}^2}{r^2+{\left({\mathrm{\σ}}_y^p\right)}^2}(y^o-{\stackrel{\‾}{y}}^p)$

所述的步骤(5)包括：

根据计算出的观测的先验集合平均和先验标准差先验集合成员以及观测误差标准差r，按照下式计算第i个集合扰动的观测增量Δy′_i:

${\mathrm{\Δy}}_i^{\′}=(\sqrt{\frac{r^2}{r^2+{\left({\mathrm{\σ}}_y^p\right)}^2}}-1)(y_i^p-{\stackrel{\‾}{y}}^p)$

所述的步骤(6)包括：

利用预先定好的集合平均的局地化因子a^EM，根据下式将集合平均的观测增量线性投影到模式状态的集合平均上：

$\mathrm{\Δ}{\stackrel{\‾}{x}}_j=\frac{{\mathrm{\ρ}}_{j,y}^{\mathrm{EM}}\·{\mathrm{cov}}_{j,y}^p}{{\left({\mathrm{\σ}}_y^p\right)}^2}\mathrm{\Δ}\stackrel{\‾}{y},$

其中，x_j代表待调整的第j个状态变量；表示观测与x_j的先验误差协方差；为集合平均的局地化函数，其计算公式为：

${\mathrm{\ρ}}_{j,y}^{\mathrm{EM}}=\mathrm{\Ω}(a^{\mathrm{EM}},b)$

其中

$\mathrm{\&Omega;}(a,b)=\left\{\begin{array}{cc}-\frac{1}{4}{\left(\frac{b}{a}\right)}^5+\frac{1}{2}{\left(\frac{b}{a}\right)}^4+\frac{5}{8}{\left(\frac{b}{a}\right)}^3-\frac{5}{3}{\left(\frac{b}{a}\right)}^2+1,& 0\&le;b\&le;a;\\ \frac{1}{12}{\left(\frac{b}{a}\right)}^5-\frac{1}{2}{\left(\frac{b}{a}\right)}^4+\frac{5}{8}{\left(\frac{b}{a}\right)}^3+\frac{5}{3}{\left(\frac{b}{a}\right)}^2-5\left(\frac{b}{a}\right)+4-\frac{2}{3}{\left(\frac{b}{a}\right)}^{-1},& a<b\&le;2a;\\ 0,& b>2a.\end{array}\right.$

其中，a为局地化因子，b为观测与待调整的状态变量之间的物理距离。

所述的步骤(7)如下：

利用预先定好的集合扰动的局地化因子a^EP，根据下式将各个集合扰动的观测增量Δy′_i线性投影到模式状态的相应的集合扰动上：

${\mathrm{\&Delta;x}}_{i,j}^{\&prime;}=\frac{{\mathrm{\&rho;}}_{j,y}^{\mathrm{EP}}\&CenterDot;{\mathrm{cov}}_{j,y}^p}{{\left({\mathrm{\&sigma;}}_y^p\right)}^2}{\mathrm{\&Delta;y}}_i^{\&prime;},$

其中为集合扰动的局地化函数，其计算公式为：

${\mathrm{\&rho;}}_{j,y}^{\mathrm{EP}}=\mathrm{\&Omega;}(a^{\mathrm{EP}},b).$

本发明的有益效果在于：

(1)对集合卡尔曼滤波中的传统局地化方法进行了改进，有效地考虑了集合平均和集合扰动所代表的不同空间尺度，显著提高了集合卡尔曼滤波的同化精度。

(2)在不引入方差膨胀的情况下，这种新的局地化方案的同化效果与引入方差膨胀的传统局地化方案得到的同化效果相当，即能够显著降低对方差膨胀的需求。

附图说明

图1为引入新的局地化方案的集合卡尔曼滤波数据同化执行流程图。

具体实施方式

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

本发明提供的是一种新的集合卡尔曼滤波局地化技术。包括对大气海洋环境实测数据进行预处理；计算观测资料处的先验观测集合成员；计算观测的先验集合平均和方差；计算集合平均的观测增量；计算各个集合扰动的观测增量；将集合平均的观测增量投影到模式状态的集合平均上；将各个集合扰动的观测增量投影到模式状态的相应集合扰动上；集合成员分析场的获得，从而更新背景场数据。本发明有效地考虑了集合平均和集合扰动所代表的不同空间尺度，显著提高了集合卡尔曼滤波的同化精度。在不引入方差膨胀的情况下，这种新的局地化方案的同化效果与引入方差膨胀的传统局地化方案得到的同化效果相当，即能够显著降低对方差膨胀的需求。本发明对实现海洋环境状态场的实时构建以及预测具有重要的实用价值。

与利用传统的集合卡尔曼滤波中的局地化技术相比，本发明的显著特征在于：对集合平均(即集合成员的平均值)和集合扰动(即集合成员与集合平均的差值)分别进行局地化，从而有效分辨集合平均和集合扰动各自代表的空间尺度。对每个观测资料，具体实现过程为：首先，利用模式的先验集合成员，观测误差标准差以及观测值，分别计算集合平均以及各个集合扰动的观测增量；其次，利用不同的局地化因子(即观测影响半径)将集合平均以及集合扰动的观测增量线性投影到周围的模式格点上，从而完成一个观测资料的同化。重复以上的步骤，从而完成一次分析步中所有观测资料的同化。此外，理论分析以及大量的同化试验表明：集合平均的局地化因子通常设置为大气或海洋的罗斯贝变形半径的特征空间尺度，而集合扰动的局地化因子应该取为介于集合平均的局地化因子和数值模式格点所能分辨的最短的波的波长之间的一个数值。本专利提出的方法可以大大地改进集合卡尔曼滤波中的传统的局地化方案的同化效果。

具体包括以下几个步骤：

步骤一：对大气海洋环境实测数据进行预处理

为了将实测环境数据同化到背景场数据库中，首先必须剔出实测数据中的无效数据和异常数据，然后将数据转化为同化所需的格式。

对每个观测资料，重复进行以下步骤：

步骤二：计算观测资料处的先验观测集合成员

对模式状态的背景集合成员，利用线性插值的方法将其投影到观测位置，从而得到观测的先验集合成员。

步骤三：计算观测的先验集合平均和方差

利用步骤二得到的观测先验集合成员，计算其先验集合平均和先验方差。

步骤四：计算集合平均的观测增量

利用计算出的观测的先验集合平均和先验方差、观测值以及观测误差标准差，计算集合平均的观测增量。

步骤五：计算各个集合扰动的观测增量

利用计算出的观测的先验集合平均和先验方差、先验集合成员以及观测误差标准差，计算各个集合扰动的观测增量。

步骤六：将集合平均的观测增量投影到模式状态的集合平均上

利用预先定好的集合平均的局地化因子，将集合平均的观测增量线性投影到模式状态的集合平均上，完成对集合平均的同化。

步骤七：将各个集合扰动的观测增量投影到模式状态的相应集合扰动上

利用预先定好的集合扰动的局地化因子，将各个集合扰动的观测增量线性投影到模式状态的相应的集合扰动上，完成对集合扰动的同化。

步骤八：集合成员分析场的获得

将由步骤七获得的集合扰动的分析场叠加到由步骤六得到集合平均的分析场上，得到最终的集合成员的分析场。

本发明提出一种新的集合卡尔曼滤波局地化同化技术，具体包括以下几个步骤：

步骤一：对实测海洋环境数据进行预处理和质量控制

为了尽可能的去除观测数据中包含的错误和不规范数据，需要首先对实测海洋环境数据进行预处理和质量控制。本专利将通过区域检验、重复深度检验、深度颠倒检验、温度数据范围检验、温度与盐度梯度检验、密度稳定性检验等六个步骤完成数据的预处理和质量控制。

对每个观测资料y^o，重复进行一下步骤：

步骤二：计算y^o的先验观测集合成员

对模式状态向量x的第i个先验集合成员，利用线性插值算子H将其投影到y^o的地理位置上，从而得到y^o的第i个先验集合成员。

步骤三：计算y^o的先验集合平均和标准差

利用步骤二得到的y^o的先验集合成员，计算其先验集合平均和先验标准差

步骤四：计算集合平均的观测增量

根据计算出的观测的先验集合平均和先验标准差观测值y^o以及观测误差标准差r，按照下式计算集合平均的观测增量

$\mathrm{\&Delta;}\stackrel{\&OverBar;}{y}=\frac{{\left({\mathrm{\&sigma;}}_y^p\right)}^2}{r^2+{\left({\mathrm{\&sigma;}}_y^p\right)}^2}(y^o-{\stackrel{\&OverBar;}{y}}^p)$

步骤五：计算各个集合扰动的观测增量Δy′_i

根据计算出的观测的先验集合平均和先验标准差先验集合成员以及观测误差标准差r，按照下式计算第i个集合扰动的观测增量Δy′_i:

${\mathrm{\&Delta;y}}_i^{\&prime;}=(\sqrt{\frac{r^2}{r^2+{\left({\mathrm{\&sigma;}}_y^p\right)}^2}}-1)(y_i^p-{\stackrel{\&OverBar;}{y}}^p)$

步骤六：将集合平均的观测增量投影到模式状态的集合平均上

利用预先定好的集合平均的局地化因子a^EM，根据下式将集合平均的观测增量线性投影到模式状态的集合平均上：

$\mathrm{\&Delta;}{\stackrel{\&OverBar;}{x}}_j=\frac{{\mathrm{\&rho;}}_{j,y}^{\mathrm{EM}}\&CenterDot;{\mathrm{cov}}_{j,y}^p}{{\left({\mathrm{\&sigma;}}_y^p\right)}^2}\mathrm{\&Delta;}\stackrel{\&OverBar;}{y},$

其中，x_j代表待调整的第j个状态变量；表示观测与x_j的先验误差协方差；为集合平均的局地化函数，其计算公式为：

${\mathrm{\&rho;}}_{j,y}^{\mathrm{EM}}=\mathrm{\&Omega;}(a^{\mathrm{EM}},b)$

其中

$\mathrm{\&Omega;}(a,b)=\left\{\begin{array}{cc}-\frac{1}{4}{\left(\frac{b}{a}\right)}^5+\frac{1}{2}{\left(\frac{b}{a}\right)}^4+\frac{5}{8}{\left(\frac{b}{a}\right)}^3-\frac{5}{3}{\left(\frac{b}{a}\right)}^2+1,& 0\&le;b\&le;a;\\ \frac{1}{12}{\left(\frac{b}{a}\right)}^5-\frac{1}{2}{\left(\frac{b}{a}\right)}^4+\frac{5}{8}{\left(\frac{b}{a}\right)}^3+\frac{5}{3}{\left(\frac{b}{a}\right)}^2-5\left(\frac{b}{a}\right)+4-\frac{2}{3}{\left(\frac{b}{a}\right)}^{-1},& a<b\&le;2a;\\ 0,& b>2a.\end{array}\right.$

其中，a为局地化因子，b为观测与待调整的状态变量之间的物理距离。

步骤七：将各个集合扰动的观测增量Δy′_i投影到模式状态的相应集合扰动上

利用预先定好的集合扰动的局地化因子a^EP，根据下式将各个集合扰动的观测增量Δy′_i线性投影到模式状态的相应的集合扰动上：

${\mathrm{\&Delta;x}}_{i,j}^{\&prime;}=\frac{{\mathrm{\&rho;}}_{j,y}^{\mathrm{EP}}\&CenterDot;{\mathrm{cov}}_{j,y}^p}{{\left({\mathrm{\&sigma;}}_y^p\right)}^2}{\mathrm{\&Delta;y}}_i^{\&prime;},$

其中为集合扰动的局地化函数，其计算公式为：

${\mathrm{\&rho;}}_{j,y}^{\mathrm{EP}}=\mathrm{\&Omega;}(a^{\mathrm{EP}},b).$

Claims (5)

1.一种集合卡尔曼滤波局地化方法，其特征在于：

(1)对大气海洋环境实测数据进行预处理

剔出实测数据中的无效数据和异常数据，然后将数据转化为同化所需的格式；

(2)对每个观测资料，计算观测资料处的先验观测集合成员：

对模式状态的背景集合成员，利用线性插值的方法投影到观测位置，得到观测的先验集合成员；

(3)计算观测的先验集合平均和方差：

利用步骤二得到的观测先验集合成员，计算其先验集合平均和先验方差；

(4)计算集合平均的观测增量：

利用计算出的观测的先验集合平均和先验方差、观测值以及观测误差标准差，计算集合平均的观测增量；

(5)计算各个集合扰动的观测增量：

利用计算出的观测的先验集合平均和先验方差、先验集合成员以及观测误差标准差，计算各个集合扰动的观测增量:；

(6)将集合平均的观测增量投影到模式状态的集合平均上：

利用:集合平均的局地化因子，将集合平均的观测增量线性投影到模式状态的集合平均上，完成对集合平均的同化；

(7)将各个集合扰动的观测增量投影到模式状态的相应集合扰动上：

利用集合扰动的局地化因子，将各个集合扰动的观测增量线性投影到模式状态的相应的集合扰动上，完成对集合扰动的同化；

(8)获得集合成员分析场

将由步骤(7)获得的集合扰动的分析场叠加到由步骤(6)得到集合平均的分析场上，得到最终的集合成员的分析场。

2.根据权利要求1所述的一种集合卡尔曼滤波局地化方法，其特征在于，所述的步骤(4)包括：

根据计算出的观测的先验集合平均和先验标准差、观测值y^o以及观测误差标准差r，按照下式计算集合平均的观测增量

$\mathrm{\&Delta;}\stackrel{\&OverBar;}{y}=\frac{{\left({\mathrm{\&sigma;}}_y^p\right)}^2}{r^2+{\left({\mathrm{\&sigma;}}_y^p\right)}^2}(y^o-{\stackrel{\&OverBar;}{y}}^p)$

3.根据权利要求1所述的一种集合卡尔曼滤波局地化方法，其特征在于，所述的步骤(5)包括：

根据计算出的观测的先验集合平均和先验标准差、先验集合成员以及观测误差标准差r，按照下式计算第i个集合扰动的观测增量Δy′_i:

$\mathrm{\&Delta;}{y}_i^{\&prime;}=(\sqrt{\frac{r^2}{r^2+{\left({\mathrm{\&sigma;}}_y^p\right)}^2}}-1)(y_i^p-{\stackrel{\&OverBar;}{y}}^p)$

4.根据权利要求1所述的一种集合卡尔曼滤波局地化方法，其特征在于，所述的步骤(6)包括：

利用预先定好的集合平均的局地化因子a^EM，根据下式将集合平均的观测增量线性投影到模式状态的集合平均上：

$\mathrm{\&Delta;}{\stackrel{\&OverBar;}{x}}_j=\frac{{\mathrm{\&rho;}}_{j,y}^{\mathrm{EM}}\&CenterDot;{\mathrm{cov}}_{j,y}^p}{{\left({\mathrm{\&sigma;}}_y^p\right)}^2}\mathrm{\&Delta;}\stackrel{\&OverBar;}{y},$

其中，x_j代表待调整的第j个状态变量；表示观测与x_j的先验误差协方差；为集合平均的局地化函数，其计算公式为：

${\mathrm{\&rho;}}_{j,y}^{\mathrm{EM}}=\mathrm{\&Omega;}(a^{\mathrm{EM}},b)$

其中

$\mathrm{\&Omega;}(a,b)=\left\{\begin{array}{cc}-\frac{1}{4}{\left(\frac{b}{a}\right)}^5+\frac{1}{2}{\left(\frac{b}{a}\right)}^4+\frac{5}{8}{\left(\frac{b}{a}\right)}^3-\frac{5}{3}{\left(\frac{b}{a}\right)}^2+1,& 0\&le;b\&le;a;\\ \frac{1}{12}{\left(\frac{b}{a}\right)}^5-\frac{1}{2}{\left(\frac{b}{a}\right)}^4+\frac{5}{8}{\left(\frac{b}{a}\right)}^3+\frac{5}{3}{\left(\frac{b}{a}\right)}^2-5\left(\frac{b}{a}\right)+4-\frac{2}{3}{\left(\frac{b}{a}\right)}^{-1},& a<b\&le;2a;\\ 0,& b>2a.\end{array}\right.$

其中，a为局地化因子，b为观测与待调整的状态变量之间的物理距离。

5.根据权利要求1所述的一种集合卡尔曼滤波局地化方法，其特征在于，所述的步骤(7)如下：

利用预先定好的集合扰动的局地化因子a^EP，根据下式将各个集合扰动的观测增量Δy′_i线性投影到模式状态的相应的集合扰动上：

$\mathrm{\&Delta;}{x}_{i,j}^{\&prime;}=\frac{{\mathrm{\&rho;}}_{j,y}^{\mathrm{EP}}\&CenterDot;{\mathrm{cov}}_{j,y}^p}{{\left({\mathrm{\&sigma;}}_y^p\right)}^2}\mathrm{\&Delta;}{y}_i^{\&prime;},$

其中为集合扰动的局地化函数，其计算公式为：

${\mathrm{\&rho;}}_{j,y}^{\mathrm{EP}}=\mathrm{\&Omega;}(a^{\mathrm{EP}},b).$