CN102004856A - 高频观测资料实时数据的快速集合卡曼滤波同化方法 - Google Patents

Abstract

本发明涉及一种高频观测资料实时数据的快速集合卡曼滤波同化方法。包括采集高频观测数据并进行质量控制；计算观测误差协方差矩阵；通过计算预报趋势，即相邻背景场的差值来获取背景场的误差协方差矩阵；利用上述协方差矩阵和背景场误差协方差矩阵，以及观测数据与海洋数值模型当前已经计算得到的背景场，对不同时刻的观测数据进行实时同化，并将经过更新的分析场赋值给下一时刻积分的初始场，继续向前预报；重复上述操作，即实现了在积分过程中对不同时刻高频观测资料的实时同化。本发明实现了高频观测资料的实时同化，提高了数据同化效率，克服了传统EnKF实施过程中需要大量集合模型同时运行的不足，规避了不收敛的问题，达到了精确数值模拟和海洋预报的目的。

Description

高频观测资料实时数据的快速集合卡曼滤波同化方法

技术领域

本发明涉及一种高频资料实时同化方法，具体涉及一种高频观测资料实时数据的快速集合卡曼滤波同化方法，属于海洋要素数值预报与海洋工程数值模拟领域。

背景技术

目前，海洋资源急需开发和管理。伴随全球气候变化而产生的海洋气候的变化，使得各种海洋灾害愈加频繁，船舶溢油事故和赤潮灾害等海洋环境污染事件也频繁发生。因此，研究全球气候变化和近海海洋灾害，对气候变化和近海海洋灾害进行预报，对海洋污染物的分布及运动路径进行预报具有重要意义。实现大范围海域精确的海洋预报，使用计算机进行海洋数值模拟是最有效的方法。然而，目前高精度、高分辨率的海洋数值模型大部分为区域海洋数值模型，这种模型严重依赖于开边界条件以及初始状态的选取，但是开边界条件和初始状态又无法客观、准确获得，只能人为调整或者给气候态平均值，适用性不佳，主观性较强，预报效果不理想。

现有的海洋预报方式主要基于数值模型。海洋数值模型是对海洋过程进行模拟计算的数学模型，在实际运用过程中都做了一定程度的简化，所以存在客观误差。误差包括模型自身的误差以及输入资料的误差。前者有结构模型、模型参数初始状态等造成的误差；而后者主要包括流速、流向等实测误差。随着积分时间的推移，预报误差会逐渐累积，从而影响预报效果。因此，在预报时通过对观测数据进行同化以对误差进行实时校正是十分必要的。另一方面，随着观测技术的发展和观测能力的提高，高精度的观测数据越来越多，比如浮标、基站、锚定、海床基、卫星和雷达等观测资料。对于流场的观测，高频地波雷达无论在空间精度还是时间精度上都有无可媲美的优势。所以，构建一个实际可行的快速实时同化方案成为必要途径，它可以充分利用高频地波雷达观测的海流数据来修正、改善海洋数值模型预报结果，有效地抑制估计误差的增长，进而提高海洋预报精度，在海洋要素数值预报与海洋工程数值模拟领域中有着重要的实用价值。对于精确海洋预报来说，这是急需解决的关键难题，但是已有的海洋数据同化方法，如卡曼滤波方法和变分方法，虽然理论方面都已经很完善，但是在实际操作过程中都有各自的局限性，无法做到高频、实时、快速同化雷达流场资料，成为精确海洋预报的瓶颈。

发明内容

本发明的目的是提供一种高频观测资料实时数据的快速集合卡曼滤波同化方法，以实现海洋数据的高频、实时、快速同化，克服现有技术的不足。

基于数值模型的海洋预报，通过对初始场的积分获得不同时刻的海洋数据。采用数据同化方法可以有效校正积分过程中的部分误差。在进行数据同化时，如果某时刻不存在观测资料，则直接进行下一步的积分运算；如果在某时刻存在相应的观测资料，则进行数据同化以校正积分过程中的误差。

本发明是在结合了加拿大快速方差法和集合卡曼滤波EnKF(Ensemble Kalman Filter)技术的基础上，基于预报误差可以被预报趋势所代替的技术构思而建立的。

高频观测资料实时数据的快速集合卡曼滤波同化方法，包括选取海洋数值模型，其特征在于还包括以下步骤：

步骤1：因高频观测资料能及时地反映海洋的状态，故首先采集高频观测资料，并通过阈值判别法和变率判别法对该观测数据进行质量控制操作，从而获得流场的观测数据矩阵D。所述的高频观测资料包括通过浮标、基站、锚定、海床基、卫星或雷达等途径获得的高频海流观测资料；所述的质量控制操作包括剔除观测中南北或者东西方向超过2米/秒的数据，剔除与平均值之差绝对值大于2倍标准差的数据，剔除前后两次观测变率大于2倍标准差的数据，剔除前后均无有效采样点的数据；若是通过地波雷达获得海流观测资料，则是同时安放两套高频地波雷达以收集雷达信号公共覆盖区域海洋表层流场观测数据；

步骤2：利用上述观测数据中观测仪器的地理位置的误差来计算观测误差协方差矩阵γγ^T，即观测数据的权重；

步骤3：在采用海洋数值模型进行预报的过程中，通过计算预报趋势，即相邻背景场的N个差值A₁′＝A₁-A₀，A₂′＝A₂-A₁，…，A_N′＝A_N-A_N-1来获取背景场的误差协方差矩阵HA′；

步骤4：利用步骤2和3中计算求得的观测误差协方差矩阵γγ^T和背景场误差协方差矩阵HA′，以及步骤1的观测数据矩阵D与海洋数值模型当前已经计算得到的背景场A_N，对不同时刻的观测数据进行实时同化，而得到经过更新的分析场A_a并将该分析场赋值给下一时刻积分的初始场，被更新后的预报状态继续向前预报；然后重复步骤1～4的操作，即实现了在积分过程中对不同时刻高频观测资料的实时同化。

为了验证步骤4中的同化结果，可以采用独立的观测数据来验证同化结果的正确性，如比较两种数据的差异，若差异变小则说明同化结果理想。

本发明所选取海洋数值模型可以是现有的任伺可以进行海流预报的海洋数值模型。

本发明改善了现有的海洋数据同化方法，提高了数据同化效率，克服了传统EnKF实施过程中需要大量集合模型同时运行的困难，且规避了传统EnKF方法在高频同化应用中出现不收敛的问题，充分利用了高频观测的表层流场资料，对整个研究海域实现了高分辨率四维、多变量数据同化，被更新的预报场不单单包含流场，还有水位、温度场以及盐度场，进而达到了进行精确数值模拟和海洋预报的目的。显然，本发明可以同化的高频观测资料可以通过浮标、基站、锚定、海床基、卫星和雷达等途径获得。

本发明结合了加拿大快速方差法和EnKF双重优点。现有的加拿大快速方差法是只通过一次预报模型的连续输出的差值来计算求得模型的背景场误差，因而最大的优点是快速、大大节省机时。EnKF技术的优点在于分析计算过程简单，不需要伴随矩阵的求解，不需要对背景场误差矩阵的反复存储和演算。但是传统EnKF在每两次数据同化间隔内，需要进行N次完整的积分运算(N为集合数)，而本发明则只需要进行1次完整的积分运算，采用初始集合最优估计做顺序积分的变化趋势，代替N个集合的依次积分做为集合扰动矩阵，在集合扰动矩阵与观测资料的误差协方差矩阵方差最小的情形下得到分析场，并以此分析场作为下一步积分的初始场。所以，本发明集合了这两种方法的优点，在大型海洋数值模型同化的实际应用中有着绝对的优势。首先，本发明实现了海洋数值模型的实时、高频同化。其次，在时机使用方面，采用了该方法的同化运算只比常规无同化的预报运算多一倍左右的运算时间，而传统意义上的EnKF同化的运算量至少是无同化运算量的N倍，因此本发明可以在保证预报精度的前提下减小计算量，提高同化系统的运算速度。而且，本发明的同化结果与实测资料的对比表明，本发明能够使温、盐、流场的结构更加趋向于观测结果，弥补了小区域模拟海流时由于驱动场、开边界和初始场精度不够造成的流场信息不准确的不足。

附图说明

图1是本发明的基本流程示意图。

其中，N为集合数，表示选取同化前N个相邻时刻的背景场，虚线为无同化的运算流程，实线为有同化的运算流程。

图2是本发明的计算背景误差协方差的流程示意图。

其中，A₀为初始场，A₁、A₂…A_N为同化前N个相邻时刻的背景场。

具体实施方式

本发明所涉及的高频观测资料包括通过浮标、基站、锚定、海床基、卫星和雷达等途径获得高频海流观测资料。

以高频地波雷达获得高频海流观测资料为例，本发明选取任意可以进行海流预报的海洋数值模型，包括以下步骤：

步骤1：同时安放两套高频地波雷达以收集雷达信号公共覆盖区域海洋表层流场观测数据，并通过阈值判别法和变率判别法对该观测数据进行质量控制操作，从而获得流场的观测数据矩阵D，即该观测数据可以用矩阵表示。高频地波雷达的资料是有误差的，这些误差或是来源于仪器本身，或是来自于安放及观测过程中，也会由于天气状况等其他因素造成。因此，观测数据是无法被直接作为真实数据来使用的，在使用前应该对观测数据进行质量控制的操作，以剔除一些不好的观测数据。所述的质量控制操作包括剔除观测中南北或者东西方向超过2米/秒的数据(属阈值判别法)，剔除与平均值之差绝对值大于2倍标准差的数据(属阈值判别法)，剔除前后两次观测变率大于2倍标准差的数据(属变率判别法)，剔除前后均无有效采样点的数据(属变率判别法)。

步骤2：利用上述观测数据中观测仪器的地理位置的误差来计算观测误差协方差矩阵γγ^T，即观测数据的权重；其中上述观测仪器的地理位置的误差为高频地波雷达自行获得的数据。现有的地波雷达测量的是海表面径向方向的流速，因此，需要同时安放两套雷达，以其信号公共覆盖区域数据求出矢量流速。由于两个雷达天线互相干扰，在计算矢量流速的过程中会出现一个与计算点的地理位置有关的误差，该误差可以设置为观测误差协方差矩阵的对角线上的值，对于非对角线上的值，可以假定观测误差在时间尺度上不相关而设置为零，或者设置为一个微小的正态分布。

步骤3：在采用海洋数值模型进行预报的过程中，通过计算预报趋势，即相邻背景场(即流场)的N个差值A₁′＝A₁-A₀，A₂′＝A₂-A₁，…，A_N′＝A_N-A_N-1来获取背景场的误差协方差矩阵HA′。在数据同化过程中，需要背景场误差协方差矩阵来计算模型预报值的权重。所谓的背景场是对该区域海洋状态的一个估计值，是通过海洋数值模型来获取的。本发明中背景场误差协方差矩阵的获取，跟加拿大快速方差法的方式致，属于现有技术。通过一个海洋数值模型连续输出获取一系列背景场，相邻背景场的差异即海洋数值模型的预报趋势，可以被看作是模型的预报误差，这样获取的误差协方差矩阵即为所需的背景场误差协方差矩阵。

步骤4：如附图2所示，利用步骤2和3中计算求得的观测误差协方差矩阵γγ^T和背景场误差协方差矩阵HA′，以及观测数据矩阵D与海洋数值模型当前已经计算得到的背景场数据A_N，对不同时刻的观测数据进行实时同化，而得到经过更新的分析场A_a并将该分析场A_a赋值给下一时刻积分的初始场，被更新后的预报状态继续向前预报；然后重复步骤1～4的操作，即实现了在积分过程中对不同时刻高频观测资料的实时同化。

具体说，如图1所示，在海洋数值模型的向前预报过程中，如果运行到有观测数据的时刻，首先根据步骤1获得观测矩阵D，然后利用步骤2和3中计算求得的观测误差协方差矩阵γγ^T和背景场误差协方差矩阵HA′，与海洋数值模型当前已经计算得到的背景场A_N，对不同时刻的观测数据进行实时同化，得到分析场A_a并更新预报状态，被更新后的预报状态作为模型的下一时刻积分的初始值继续向前预报，直到再次遇到观测数据，重复以上步骤1～4的操作。

为了验证步骤4中的同化结果，可采用独立的观测数据来验证同化结果的正确性，如比较两种数据的差异，若差异变小则说明同化结果理想，以确保下一次积分的准确性。

上述步骤2和3并无先后关系，在进行数据处理时可以同时操作以提高效率。本发明的计算公式如下：

A_a＝A+A′A′^TH^T(HA′A′^TH^T+γγ^T)^-1D′，

其中，A′∈R^n×N是由相邻背景场的N个差值构成的集合扰动矩阵；A∈R^n×N是集合扰动矩阵A′与取平均后得到的

相减，再加上海洋数值模型运算的背景场A_N；n为向量A_N的维数，即向量A_N含有n个网格点数据；m为观测值的数目，H为将海洋数值模型运算得到的n个网格点数据投影至m个观测点位置的观测算子，不需要单独计算；D∈R^m×N为观测数据矩阵，D′＝D-HA∈R^m×N；γ∈R^n×N为观测扰动值集合，γγ^T为观测误差协方差矩阵，N为集合数，即选取相邻时刻背景场的个数。

本发明中，最核心的问题是简化N次集合预报的过程，将同化时间步长内N次运算结果分别与前一次运算结果相减的差值矩阵作为集合扰动矩阵，以此计算分析场A_a。一个同化时间步长内的具体计算过程如下：

①A₀→A₁→A₂→…→A_N-1→A_N∈R^n×1

②A₁′＝A₁-A₀，A₂′＝A₂-A₁，…，A_N′＝A_N-A_N-1∈R^n×1

A′＝(A₁′，A₂′，…，A_N′)∈R^n×N， $\stackrel{\‾}{A^{\′}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\left({A_i}^{\′}\right)\∈R^{n\×1}$

③ $B^{\′}=A^{\′}-\stackrel{\‾}{A^{\′}}\∈R^{n\×N}$

④A＝B′+A_N∈R^n×N

⑤ ${\mathrm{HA}}^{\′}=\mathrm{HA}-\stackrel{\‾}{\mathrm{HA}}\∈R^{m\×N}$

A₀为初始场，A₁、A₂…A_N为同化前N个相邻时刻的背景场，自然地，N即相邻时刻背景场的个数，也就是集合数；如图2所示，相邻的两个背景场相减得到集合扰动矩阵A′。取平均得到

后，集合扰动矩阵A′与相减得到B′。B′与该时刻同化前海洋数值模型的运算结果A_N的加和作为集合矩阵A。HA为集合中向量在观测点的投影，每一列集合矩阵A投影到观测点位置均得到的m个值。最终根据根据上述的计算结果HA′、γγ^T、D与已知的A_N、采用Evensen提出的EnKF方法计算分析场A_a。

Claims (8)

1.高频观测资料实时数据的快速集合卡曼滤波同化方法，包括选取海洋数值模型，其特征在于还包括以下步骤：

步骤1：首先采集海流的高频观测资料，并通过阈值判别法和变率判别法对该观测数据进行质量控制操作，从而获得流场的观测数据矩阵D；

步骤2：利用上述观测数据中观测仪器的地理位置的误差来获得观测误差协方差矩阵γγ^T，即观测数据的权重；

步骤3：在采用上述海洋数值模型进行预报的过程中，通过计算预报趋势，即相邻背景场的N个差值A₁′＝A₁-A₀，A₂′＝A₂-A₁，…，A_N′＝A_N-A_N-1来获取背景场的误差协方差矩阵HA′；

步骤4：利用步骤2和3中计算求得的观测误差协方差矩阵γγ^T和背景场误差协方差矩阵HA′，以及步骤1的观测数据矩阵D与海洋数值模型当前已经计算得到的背景场A_N，对不同时刻的观测数据进行实时同化，而得到经过更新的分析场A_a并将该分析场A_a赋值给下一时刻积分的初始场，被更新后的预报状态继续向前预报；然后重复步骤1～4的操作，即实现了在积分过程中对不同时刻高频观测资料的实时同化。

2.如权利要求1所述的同化方法，其特征在于上述步骤1中的海流的高频观测资料包括通过浮标、基站、锚定、海床基、卫星或雷达获得的海流高频观测资料。

3.如权利要求2所述的同化方法，其特征在于上述通过雷达获得海流的高频观测资料是同时安放两套高频地波雷达以收集雷达信号公共覆盖区域海洋表层流场观测数据。

4.如权利要求1所述的同化方法，其特征在于上述步骤1中所述的通过阈值判别法和变率判别法对该观测数据进行质量控制操作，包括剔除观测中南北或者东西方向超过2米/秒的数据，剔除与平均值之差绝对值大于2倍标准差的数据，剔除前后两次观测变率大于2倍标准差的数据，剔除前后均无有效采样点的数据。

5.如权利要求1所述的同化方法，其特征在于上述步骤3中背景场误差协方差矩阵HA′是以加拿大快速方差法获取。

6.如权利要求1所述的同化方法，其特征在于在上述步骤4中对不同时刻的观测数据进行实时同化而得到经过更新的分析场，是采用下列计算公式进行：

A_a＝A+A′A′^TH^T(HA′A′^TH^T+γγ^T)^-1D′，

其中，其中，A′∈R^n×N是由相邻背景场的N个差值构成的集合扰动矩阵；A∈R^n×N是集合扰动矩阵A′与取平均后得到的

相减，再加上海洋数值模型运算的背景场A_N；n为向量A_N的维数，即向量A_N含有n个网格点数据；m为观测值的数目，H为将海洋数值模型运算得到的n个网格点数据投影至m个观测点位置的观测算子，不需要单独计算；D∈R^m×N为观测数据矩阵，D′＝D-HA∈R^m×N。

7.如权利要求1所述的同化方法，其特征在于在上述步骤4的对不同时刻的观测数据进行实时同化之后，采用独立的观测数据来验证同化结果的正确性，比较两种数据的差异，若差异变小则说明同化结果理想。

8.如权利要求1所述的同化方法，其特征在于所选取的海洋数值模型是现有的任何可以进行海流预报的海洋数值模型。

Images