CN104899448A - 一种集合卡尔曼滤波静态局地化方案的自适应补偿方法 - Google Patents

Abstract

本发明属于实测海洋环境数据同化领域，具体涉及一种集合卡尔曼滤波静态局地化方案的自适应补偿方法。本发明包括：对大气海洋环境实测数据进行预处理；对每个观测资料，根据实时的观测系统，利用不同观测要素的观测误差、观测数目以及显著性水平，计算后续需要使用的阈值；顺序地同化所有的观测资料；计算其集合平均和集合扰动；计算观测余量并更新集合成员。本发明对集合卡尔曼滤波中的静态局地化方法进行了改进，有效地提取了观测信息中静态局地化方法未能提取的多尺度信息，显著提高了集合卡尔曼滤波的同化精度。

Description

一种集合卡尔曼滤波静态局地化方案的自适应补偿方法

技术领域

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

背景技术

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

本发明主要是对集合卡尔曼滤波中传统的静态局地化方案进行改进，提高集合卡尔曼滤波的同化效果。根据前述，有限的集合样本会给集合卡尔曼滤波数据同化引入显著的抽样误差，比如模式状态先验方差的低估以及模式状态与远距离观测之间的伪相关性等。解决这个问题的一个有效途径是在同化过程中引入局地化方法，即一个观测仅允许影响周围一定范围内的模式格点，且距离越远影响越小。传统的集合卡尔曼滤波局地化方案是一种静态方法，即预先给定一个观测影响半径(即局地化因子)，在整个同化过程中均采用这个局地化因子。大量的文献记载以及试验证明，集合卡尔曼滤波的同化效果对局地化因子的选取非常敏感。这种静态的局地化方案存在明显的局限性，因为对于实际的三维海洋数值模式来说，要找到一个最优的局地化因子是非常耗时的，几乎是不可能的。因此，针对静态局地化方法存在的弊端，引入一种补偿方法，从而在保证同化精度的同时降低同化效果对局地化因子的敏感性是非常有意义和有实用价值的。

发明内容

本发明的目的在于提供一种提高精度的集合卡尔曼滤波静态局地化方案的自适应补偿方法。

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

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

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

(2)对每个观测资料，根据实时的观测系统，利用不同观测要素的观测误差、观测数目以及显著性水平，计算后续需要使用的阈值；

(3)顺序地同化所有的观测资料：

利用静态局地化因子以及标准的集合卡尔曼滤波方法，对每个观测资料进行同化，对模式状态的先验集合成员进行更新，得到集合卡尔曼滤波的后验集合成员；

(4)利用步骤(3)中的集合卡尔曼滤波的后验集合成员，计算其集合平均和集合扰动；

(5)计算观测余量并更新集合成员：

将步骤(4)得到的后验集合平均投影到观测位置，从而得到后验观测值，从原始观测值中减去后验观测值，得到观测余量，并计算原始观测值和后验观测值之间的均方根误差RMSE，如果RMSE大于步骤(2)得到的阈值，则：

(5.1)提取观测余量中的多尺度信息：

利用多重网格分析技术(MGA)提取观测余量中的多尺度信息，从而得到MGA的分析场；

(5.2)对后验集合平均进行更新：

将MGA得到的分析场叠加到步骤四得到集合卡尔曼滤波的后验集合平均上，得到最终的集合平均的分析场；

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

所述的步骤(2)包括：

根据特定的观测系统，按照下式计算自适应补偿机制所需的阈值(θ)：

$\mathrm{\θ}(\mathrm{\α},K,r)\≡\sqrt{\frac{{\mathrm{\χ}}_{1-\mathrm{\α}}^2\left(K\right)}{K}}\·r,$

其中，α为显著性水平(通常取为0.01)，K为观测数，r为观测误差标准差，为自由度为K的卡方分布的(1-α)上侧分位数。

所述的步骤(3)包括：

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

(3.1)计算y^o的先验观测集合成员

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

(3.2)计算y^o的观测增量

根据第(3.1)步计算的y^o的先验集合、观测值(y^o)以及观测误差标准差(r)，按照下式计算y^o的观测增量Δy_i：

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

其中，表示y^o的第i个先验集合成员，表示y^o的先验标准差，表示y^o的先验集合平均；

(3.3)对状态变量进行更新

基于一个静态局地化函数ρ_j,y，根据计算出的y^o的观测增量，按照下式对状态变量进行更新：

${\mathrm{\Δx}}_{i,j}=\frac{{\mathrm{\ρ}}_{j,y}{\mathrm{cov}}_{j,y}^p}{{\left({\mathrm{\σ}}_y^p\right)}^2}{\mathrm{\Δy}}_i,$

其中j表示第j个状态变量，表示观测与x_j的先验误差协方差；ρ_j,y的计算公式为：

${\mathrm{\&rho;}}_{j,y}=\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^o与状态变量x_j之间的物理距离。

本发明的有益效果在于：

(1)对集合卡尔曼滤波中的静态局地化方法进行了改进，有效地提取了观测信息中静态局地化方法未能提取的多尺度信息，显著提高了集合卡尔曼滤波的同化精度。

(2)降低了集合卡尔曼滤波对静态局地化因子的强依赖性。

附图说明

图1一种集合卡尔曼滤波静态局地化方案的自适应补偿方法的执行流程图。

具体实施方式

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

本发明提供的是一种集合卡尔曼滤波静态局地化方案的自适应补偿机制。包括对大气海洋环境实测数据进行预处理；根据实时的观测系统，利用不同观测要素的观测误差、观测数目以及显著性水平，计算后续需要使用的阈值；利用集合卡尔曼滤波方法顺序地同化所有的观测资料；计算集合卡尔曼滤波的后验集合平均和集合扰动；计算观测余量并更新集合成员。本发明对集合卡尔曼滤波中的静态局地化方法进行了改进，有效地提取了观测信息中静态局地化方法未能提取的多尺度信息，显著提高了集合卡尔曼滤波的同化精度。降低了集合卡尔曼滤波对静态局地化因子的强依赖性。本发明对实现海洋环境状态场的实时构建以及预测具有重要的实用价值。

本发明针对集合卡尔曼滤波中的静态局地化方案提出一种自适应补偿方法。与传统的静态局地化技术相比，本发明的显著：该方法能够自适应地提取观测资料中静态局地化方案未能提取的观测信息。对每个分析步，具体实现过程为：首先，基于一个静态局地化因子，利用观测资料对模式的先验集合成员进行调整；其次，将模式状态的后验集合平均值投影到观测位置上，得到观测的后验值；第三，从原始观测值中减去上述后验观测值，得到观测余量，并计算原始观测与后验观测之间的均方根误差；第四，当上述均方根误差大于某个阈值时，则利用一种多重网格技术(MGA)提取观测余量中的多尺度信息，从而得到相应的分析场；第五，当上述均方根误差大于某个阈值时，将MGA的分析场叠加到第一步得到后验集合平均上得到最终的集合卡尔曼滤波的集合平均分析场。本专利提出的方法可以大大地改进集合卡尔曼滤波中的传统的静态局地化方案的同化效果。

每个分析步的同化过程具体包括以下几个步骤：

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

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

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

步骤二：根据实时的观测系统，利用不同观测要素的观测误差、观测数目以及显著性水平，计算后续需要使用的阈值。

步骤三：顺序地同化所有的观测资料

利用一个既定的静态局地化因子以及标准的集合卡尔曼滤波方法，对每个观测资料进行同化，从而对模式状态的先验集合成员进行更新，得到集合卡尔曼滤波的后验集合成员。

步骤四：利用步骤三中的集合卡尔曼滤波的后验集合成员，计算其集合平均和集合扰动。

步骤五：计算观测余量并更新集合成员

将步骤四得到的后验集合平均投影到观测位置，从而得到后验观测值。从原始观测值中减去后验观测值，得到观测余量，并计算原始观测值和后验观测值之间的均方根误差(RMSE)。如果RMSE大于步骤二得到的阈值，则：

步骤5.1：提取观测余量中的多尺度信息

利用多重网格分析技术(MGA)提取观测余量中的多尺度信息，从而得到MGA的分析场。

步骤5.2：对后验集合平均进行更新

将MGA得到的分析场叠加到步骤四得到集合卡尔曼滤波的后验集合平均上，得到最终的集合平均的分析场。

步骤5.3：将最终的集合平均分析场叠加到步骤四得到的集合扰动上，得到最终的集合成员的分析场。

本发明提出一种集合卡尔曼滤波静态局地化方案的自适应补偿方法，具体包括以下几个步骤：

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

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

步骤二：计算自适应补偿方法所需的阈值

根据特定的观测系统，按照下式计算自适应补偿方法所需的阈值(θ)：

$\mathrm{\&theta;}(\mathrm{\&alpha;},K,r)\&equiv;\sqrt{\frac{{\mathrm{\&chi;}}_{1-\mathrm{\&alpha;}}^2\left(K\right)}{K}}\&CenterDot;r,$

其中，α为显著性水平(通常取为0.01)，K为观测数，r为观测误差标准差，为自由度为K的卡方分布的(1-α)上侧分位数。

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

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

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

步骤3.2：计算y^o的观测增量

根据第3.1步计算的y^o的先验集合、观测值(y^o)以及观测误差标准差(r)，按照下式计算y^o的观测增量Δy_i:

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

其中，表示y^o的第i个先验集合成员，表示y^o的先验标准差，表示y^o的先验集合平均。

步骤3.3：对状态变量进行更新

基于一个静态局地化函数ρ_j,y，根据计算出的y^o的观测增量，按照下式对状态变量进行更新:

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

其中j表示第j个状态变量，表示观测与x_j的先验误差协方差；ρ_j,y的计算公式为：

${\mathrm{\&rho;}}_{j,y}=\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^o与状态变量x_j之间的物理距离。

步骤四：计算后验观测值

计算步骤三得到的模式状态的后验集合平均并利用观测算子H将其投影到观测位置，得到后验观测值

步骤五：计算观测余量

从原始观测向量y^o中减去后验观测值从而得到观测余量此外，计算原始观测与后验观测之间的均方根误差RMSE。当RMSE大于θ时，则顺序进行下列操作：

步骤5.1：利用一种多重网格分析技术(MGA)，提取观测余量y^res中的多尺度信息，得到MGA的分析场；

步骤5.2：将MGA的分析场叠加到步骤三得到的集合平均的分析场上，得到集合平均最终的分析场；

步骤5.3：将步骤三得到的集合扰动的分析场叠加到步骤5.2得到的集合平均的最终分析场上，得到最终的集合分析场。

Claims (3)

1.一种集合卡尔曼滤波静态局地化方案的自适应补偿方法，其特征是：

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

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

(2)对每个观测资料，根据实时的观测系统，利用不同观测要素的观测误差、观测数目以及显著性水平，计算后续需要使用的阈值；

(3)顺序地同化所有的观测资料：

利用静态局地化因子以及标准的集合卡尔曼滤波方法，对每个观测资料进行同化，对模式状态的先验集合成员进行更新，得到集合卡尔曼滤波的后验集合成员；

(4)利用步骤(3)中的集合卡尔曼滤波的后验集合成员，计算其集合平均和集合扰动；

(5)计算观测余量并更新集合成员：

将步骤(4)得到的后验集合平均投影到观测位置，从而得到后验观测值，从原始观测值中减去后验观测值，得到观测余量，并计算原始观测值和后验观测值之间的均方根误差RMSE，如果RMSE大于步骤(2)得到的阈值，则：

(5.1)提取观测余量中的多尺度信息：

利用多重网格分析技术(MGA)提取观测余量中的多尺度信息，从而得到MGA的分析场；

(5.2)对后验集合平均进行更新：

将MGA得到的分析场叠加到步骤四得到集合卡尔曼滤波的后验集合平均上，得到最终的集合平均的分析场；

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

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

根据特定的观测系统，按照下式计算自适应补偿机制所需的阈值(θ)：

$\mathrm{\&theta;}(\mathrm{\&alpha;},K,r)\&equiv;\sqrt{\frac{{\mathrm{\&chi;}}_{1-\mathrm{\&alpha;}}^2\left(K\right)}{K}}\&CenterDot;r,$

其中，α为显著性水平(通常取为0.01)，K为观测数，r为观测误差标准差，为自由度为K的卡方分布的(1-α)上侧分位数。

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

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

(3.1)计算y^o的先验观测集合成员

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

(3.2)计算y^o的观测增量

根据第(3.1)步计算的y^o的先验集合、观测值(y^o)以及观测误差标准差(r)，按照下式计算y^o的观测增量Δy_i：

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

其中，表示y^o的第i个先验集合成员，表示y^o的先验标准差，表示y^o的先验集合平均；

(3.3)对状态变量进行更新

基于一个静态局地化函数ρ_j,y，根据计算出的y^o的观测增量，按照下式对状态变量进行更新：

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

其中j表示第j个状态变量，表示观测与x_j的先验误差协方差；ρ_j,y的计算公式为：

${\mathrm{\&rho;}}_{j,y}=\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^o与状态变量x_j之间的物理距离。