一种基于集合卡尔曼滤波的潮汐数值模型水深估计方法

1.本发明涉及海洋数据同化技术，特别涉及一种基于集合卡尔曼滤波(ensemble adjustment kalman filter，eakf)的潮汐数值模型水深估计方法，主要应用于海洋潮汐数值模拟和预报。


背景技术：

2.随着计算机技术的不断发展，海洋数值模型越来越成为人们研究和预测海洋的重要工具。对于所有的海洋数值模式而言，参数化方案以及精确的参数值的给定对数值模拟结果有着重要的影响。在海洋数值模式调试的过程中，参数值的给定通常采用试错的方法，以使模拟结果接近观测。由于海洋数值模型的复杂性，这是一个极其耗费计算资源与人力资源的过程。因此，替代这种主观调试方法的途径之一是采用海洋数据同化技术，通过将观测数据同化到海洋数值模型中，在模型状态场调整的同时，对参数进行估计，从而得到一个合理的参数值。这样做的同时，也可以在一定程度上得以缓解试错法在海洋数值模型调试中所带来的繁重任务。利用基于伴随模型的四维变分同化方法与集合卡尔曼滤波方法均可以开展海洋数值模型中参数的最优估计。
3.以潮汐潮流的数值模拟为例，模型中不确定的参数包括水深、底摩擦系数和开边界条件等。此前，研究者们大多采用伴随方法来进行潮汐潮流模拟中这些参数的优化，从早期基于线性浅水方程模型的理想数据同化试验，到后来基于非线性二维和三维潮波数值模型的实际数据同化试验，均取得了很好的研究成果。然而，伴随方法有其局限性：对于不同的模式，伴随方法需要编写对应的伴随模式来求取模式的目标函数梯度，因此该方法工作量高、可移植性差。相比较而言，集合卡尔曼滤波方法不需要编写伴随模式，可实现并行化计算，工作量低、可移植性强。然而，基于集合卡尔曼滤波方法开展潮汐潮流模拟中的参数估计研究尚较少。目前多数研究着眼于近岸、河流和潮汐汊道的水深参数最优估计。而对于陆架浅海的水深估计研究尚未见相关报道。


技术实现要素：

4.本发明针对陆架浅海潮汐数值模型，主要目的是通过水深参数估计方案设计，实现eakf方法在浅海潮汐数值模型中的水深估计，进而为潮汐潮流的分析和预报提供一种新的参数优化方法，为提高潮汐潮流数值模拟精度提供技术支持，预期在潮汐潮流数值预报中具有良好的应用前景。
5.本发明所采用的技术方案是：一种基于集合卡尔曼滤波的潮汐数值模型水深估计方法，包括以下步骤：
6.步骤1，海洋数值模式水深参数的确定：根据海区海底地形和模式的敏感性分析结果，以水深大小对海区进行分区，进行敏感性分析，确定参数化方案；
7.步骤2，生成集合，进行数值模拟：在各个分区上分别设置一水深增量，将各个水深增量作为参数；通过在各个水深增量参数上叠加无偏的高斯随机数生成水深增量参数集
合，再根据水深增量参数集合生成水深参数集合，将水深参数集合代入海洋数值模式中进行自由积分直到海洋数值模式稳定；
8.步骤3，采用强化的参数校正数据同化方法(data assimilation scheme for enhancive parameter correction，daepc)进行同化：先进行海洋数值模式状态估计，在海洋数值模式状态达到准平衡状态后，再启动对海洋数值模式参数的估计，其中，采用集合调整卡尔曼滤波方法进行海洋数值模式状态估计和海洋数值模式参数估计；
9.步骤4，同化结束后，使海洋数值模式以同化后的参数自由积分至稳定，以获取优化后的海洋数值模式状态变量；
10.步骤5，将海洋数值模式状态变量进行调和分析，得到优化后的潮汐调和常数分析结果，用于潮汐预报。
11.进一步地，步骤2中，对于海洋数值模式中的任一网格点，该网格点中的水深增量参数集合表示为{
△
x1,
△
x2,
…
,
△
x
i
,
…
,
△
x
n
}，i＝1,2,
…
,n，其中，
△
x
i
表示该网格点处增量参数的第i个集合成员，n表示该网格点中水深增量参数的集合成员数量；
12.根据水深增量参数集合生成水深参数集合{ζ1,ζ2,
…
,ζ
i
,
…
,ζ
n
}，其中，ζ
i
根据公式(1)计算得到：
13.ζ
i
＝ζ
data
+
△
x
i
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
14.式中，ζ
i
表示在该网格点处的第i个水深集合成员，ζ
data
为该网格点的地形数据。
15.本发明的有益效果是：本发明一种基于集合卡尔曼滤波的潮汐数值模型水深估计方法，一方面，eakf作为一种顺序同化方法，不需要编写伴随模式和切线性模式，相比于伴随方法实现难度更小，且能有效地保留非线性信息，并易于实现并行计算，进而减少运算时长。另一方面，该方法可为提高陆架浅海区域潮汐分析和预报的精度提供新的思路。
附图说明
16.图1为本发明基于集合卡尔曼滤波的潮汐数值模型水深估计方法流程图。
17.图2为基于有偏水深数据获得的黄渤海m2分潮水位时间平均均方根误差(单位：m)。
18.图3为同化试验中黄渤海m2分潮水位均方根误差(单位：m)时间序列。
19.图4为利用优化后的水深参数获得的黄渤海m2分潮同潮图(虚线为振幅，单位：cm；实线为迟角，单位：
°
)。
20.其中，所有附图中的数值模拟结果均基于广义坐标系统的美国普林斯顿大学海洋模式(princeton ocean model with generalized coordinate system，简称pomgcs)。
具体实施方式
21.为能进一步了解本发明的发明内容、特点及功效，兹例举以下实施例，并配合附图详细说明如下：
22.如附图1所示，一种基于集合卡尔曼滤波的潮汐数值模型水深估计方法，为一种根据海区地形、设计水深的参数化方案，通过eakf将观测数据同化到潮汐数值模型中，对水深参数进行最优估计，提高潮汐模拟精度的方法。包括以下步骤：
23.(1)海洋数值模式水深参数估计方案的确定
24.根据海区海底地形和模式的敏感性分析结果，以水深大小划分海区，进行敏感性分析，确定参数化方案。
25.基于图2中绘制的黄渤海m2分潮的敏感性分析结果(用水位数值模拟结果的时间平均均方根误差表示)，对确定模式水深参数的流程进行详细描述。首先，根据所选海区的海底地形，给定一参数化方案：图2中，将海区按照深度分为0
‑
20m、20
‑
40m、40
‑
60m、>60m四个区域。根据水深数据的精度，在上述四个区域的水深数据基础上，分别加0.5m、2m、4m和6m，生成相应的有偏水深数据，并保持模式其它参数(模式里除了水深参数之外的其余全部参数)不变，进行10天自由积分。在最后一个周期里，将水位计算结果与原模式结果进行对比，计算各网格点水位的时间平均均方根误差。通过敏感性分析，能够看出参数误差可以引起足够的模式误差，说明该参数化方案满足数据同化实验的要求。
26.(2)集合的生成
27.在各个分区上分别设置一水深增量，将各个水深增量作为参数；通过在各个水深增量参数上叠加无偏的高斯随机数生成水深增量参数集合，再根据水深增量参数集合生成水深参数集合，将水深参数集合代入海洋数值模式中进行自由积分直到海洋数值模式稳定。
28.对于海洋数值模式中的任一网格点，该网格点中的水深增量参数集合表示为{
△
x1,
△
x2,
…
,
△
x
i
,
…
,
△
x
n
}，i＝1,2,
…
,n，其中，
△
x
i
表示该网格点处增量参数的第i个集合成员，n表示该网格点中水深增量参数的集合成员数量；
29.根据水深增量参数集合生成水深参数集合{ζ1,ζ2,
…
,ζ
i
,
…
,ζ
n
}，其中，ζ
i
根据公式(1)计算得到：
30.ζ
i
＝ζ
data
+
△
x
i
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
31.式中，ζ
i
表示在该网格点处的第i个水深集合成员，ζ
data
为该网格点的地形数据。
32.(3)同化试验
33.采用强化的参数校正数据同化方法进行同化：该方法的实现方式是在海洋数值模式状态场的估计达到准平衡态后，再启动对海洋数值模式参数的估计。这样可以先减小海洋数值模式的状态误差，使海洋数值模式误差以参数误差为主导，进而强化参数调整过程中的信噪比，从而提高参数估计的效果。对于本发明，海洋数值模式状态包括水位和潮流，海洋数值模式参数为水深。采用eakf方法进行状态和参数估计，根据daepc方法的实现方式，先进行水位和潮流的状态估计，当水位和潮流达到准平衡状态时，再加入水深的参数估计。
34.基于最小二乘框架的eakf方法的原理如下所示：
35.以某一观测y
o
为例，该方法的实施包括如下两个步骤：
36.第一步，计算观测增量：
[0037][0038]
式中，为第i个集合成员在观测点上的观测增量；σ
o
为观测误差的标准差；为集合在该观测点上的先验(同化前)误差标准差；为第i个集合成员在该观测点上的先
验值；为投影到该观测点上的集合平均值。
[0039]
第二步，将观测增量投影到模式网格点上，获得集合成员的更新值：
[0040][0041]
式中，为第i个先验集合成员；为第i个后验集合成员；为状态集合与观测之间的协方差；r
x,y
为局地化因子，采用gaspari
‑
cohn函数形式：
[0042][0043]
式中，a为经验给定的局地化半径；b为观测点与状态变量格点之间的距离。
[0044]
为保证集合的离散度，避免滤波发散。对于状态估计，本发明引入静态膨胀方案。即确定一个常数膨胀因子，对各个集合成员相对于集合平均的扰动进行膨胀，用以调整集合的离散度，从而避免集合离散度降低。而对于参数估计，则引入条件静态乘法膨胀方案，即通过判断当前时刻集合方差的大小，决定是否对其进行参数膨胀。条件静态乘法膨胀公式如下：
[0045][0046]
式中，β为膨胀前的参数集合；为膨胀后的参数集合；为参数的集合平均；α0为根据经验设定的参数膨胀系数；σ0为初始集合的标准差；σ
t
为t时刻集合的标准差。
[0047]
基于上述方法，进行数据同化。在同化的实现过程中，采用daepc方法：首先仅进行状态估计至准平衡状态。待海洋数值模式稳定，加入参数估计并同化至海洋数值模式再次稳定。
[0048]
(4)获取优化的调和常数
[0049]
同化结束后，使海洋数值模式以同化后的参数自由积分至稳定，以获取优化后的海洋数值模式状态变量。将海洋数值模式状态变量进行调和分析，得到优化后的潮汐调和常数分析结果，可以用于潮汐预报。本实施例中，将水位进行调和分析，获取m2分潮的调和常数。
[0050]
图3为同化试验中黄渤海m2分潮水位均方根误差时间序列。图4为利用优化后的水深参数获得的黄渤海m2分潮同潮图。
[0051]
尽管上面结合附图对本发明的优选实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，并不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可以做出很多形式，这些均属于本发明的保护范围之内。