一种基于数据同化的近岸水质数据再分析方法与流程

1.本发明创造属于水环境数据再分析和预报技术领域，尤其是涉及一种基于数据同化的近岸水质数据再分析方法。


背景技术：

2.随着计算机技术的发展，数值模拟已经成为了水动力、水环境分析和预报的重要手段。市面上已有delft3d，swan等开源软件和mike21、efdc等商业软件。但仅依靠数值模型对较复杂的水动力过程进行分析和预报往往不能得到理想的效果：一方面模型的离散格式和模型方程的近似带来了一定的模型误差；另一方面人为给定的底摩擦系数、地形和边界条件等模型参数也对模型的计算精度产生了较大影响。水环境数值模拟中，目前广泛应用的数据同化方法主要包括最优插值、四维变分、粒子滤波和集合卡尔曼滤波（ensemble kalman filter，enkf）等。其中，四维变分和enkf能够在调整模型状态场的同时对模型参数进行估计，从而获得更为准确的参数值。其中，enkf方法能够量化模型的不确定性和误差，并随着时间步的推进不断对模式进行优化，相比于同样是目前研究热点的四维变分方法，enkf方法实现相对简单，且更加适用于强非线性的数学模型。目前已有采用enkf方法开展水环境数值模拟的相关研究，然而，多数研究中采用的是传统的enkf方法，需要的集合成员数较大，模型运算时间较长，除此之外，由于enkf方法在实现过程中需要给定一观测扰动集合，在实际应用中，生成该集合时会产生一定的抽样误差，从而导致同化精度的降低，因此，有必要进行优化改进。


技术实现要素：

3.有鉴于此，本发明创造旨在克服现有技术中的缺陷，提出一种基于数据同化的近岸水质数据再分析方法，实现集合调整卡尔曼滤波（ensemble adjustment kalman filter，eakf）方法对水环境模拟中水质的状态估计，进而为水环境的分析和预报提供一种新的优化方法，为提高水环境数值再分析和预报精度提供技术支持。
4.为达到上述目的，本发明创造的技术方案是这样实现的：一种基于数据同化的近岸水质数据再分析方法，包括如下步骤：基于文献资料和实测数据进行水动力模型的设置：水动力模型可采用冷启动，即初始条件可全部设置为0，也可基于已有水位、流速数据进行设置；可采用验潮站观测数据、卫星数据或再分析数据集等给定开边界处的水位和流速，波浪开边界采用侧向边界条件；由于研究区域范围较小，可假设区域内底摩擦处处相同，即曼宁系数可设为常数；模型水深可采用卫星测深数据，根据模型分辨率进行插值生成；风场可采用ecwmf等再分析数据。通过在模型原有曼宁系数的基础上叠加无偏的高斯随机数生成参数集合，使模型自由积分至稳定；采用强化的参数校正数据同化方法进行数据同化，该数据同化方法包括：先进行潮位、潮流和波高的状态估计，待模型达到准平衡状态后，同时进行状态估计和对曼宁系数
的参数估计；反复调试数据同化算法中的局地化半径和膨胀系数，以获取同化效果最优的模型状态和参数，该过程通过人为设置“真实”解并基于“真实”解生成观测场，测试不同局地化半径和膨胀系数下同化结果与“真实”解的误差，最终确定一组最优的同化参数；水环境模型的初始条件由对应时间步的水动力模型的数值模拟结果进行插值生成，水环境模型的动力边界条件采用水动力模型的潮位、潮流模拟结果，水质边界条件采用零梯度边界条件，通过在水质扩散系数的基础上叠加无偏的高斯随机数生成参数集合，水质指标包括但不限于溶解氧、化学需氧量、生物需氧量、氨氮、总磷、总氮、硝酸盐氮、叶绿素等；基于上述生成的水环境模型边界条件和实测温度、盐度等模型参数条件，进行水环境数值模拟，将曼宁系数设置为水动力模型中优化后的结果，进行模型自由积分；通过与水动力模型形式相同的理想试验测试模型的最优同化参数，再采用最优同化参数，基于eakf方法对水动力和水质浓度同时进行状态估计并针对水质扩散系数开展参数估计，获取优化后的水动力和水质数据。
5.进一步，水动力模型、水环境模型的风场均采用实测或再分析资料；水环境模型的水质浓度选择实测数据。
6.进一步，水环境模型参数中曼宁系数采用eakf优化后的结果。
7.进一步，水动力模型的运算中，将模型设置为冷启动，给定各时刻潮位、潮流数据作为边界条件，根据海区具体海底地貌，曼宁系数取值不同，可参照相关文献或已有研究结果设定模型的曼宁系数值。
8.进一步，水环境模型的运算中，将模型设置为热启动。
9.进一步，实测数据包括地形、风场、水质、海区温度、海区盐度等初始场和数据同化观测数据，以及潮位、潮流边界数据。
10.进一步，文献资料包括ecwmf的风场资料，以及验潮站的潮位、潮流数据，以及nao.99b的潮位、潮流数据。
11.进一步，水动力模型的边界条件直接采用潮位站数据。
12.相对于现有技术，本发明创造具有以下优势：本发明方法将实测数据和数值模拟相结合，通过基于最小二乘的eakf方法对模型参数和状态变量进行订正，最终提高模型模拟精度，获得较高精度水环境数据的方法，本方法主要应用于水环境数据的再分析和预报，基于观测数据，通过数据同化，将观测信息与数值模型相结合，以此优化水环境模型的模拟结果，从而提供精度较高的水环境再分析和预报数据。
附图说明
13.构成本发明创造的一部分的附图用来提供对本发明创造的进一步理解，本发明创造的示意性实施例及其说明用于解释本发明创造，并不构成对本发明创造的不当限定。在附图中：图1为本发明基于数据同化的近岸水质数据再分析方法的示意图；图2为本发明中基于集合调整卡尔曼滤波的水环境数据再分析和预报方法的流程
图。
具体实施方式
14.需要说明的是，在不冲突的情况下，本发明创造中的实施例及实施例中的特征可以相互组合。
15.在本发明创造的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明创造和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明创造的限制。此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”等的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明创造的描述中，除非另有说明，“多个”的含义是两个或两个以上。
16.在本发明创造的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以通过具体情况理解上述术语在本发明创造中的具体含义。
17.下面将参考附图并结合实施例来详细说明本发明创造。
18.一种基于数据同化的近岸水质数据再分析方法，如图1和图2所示，包括如下步骤：以氨氮数据再分析过程为例，说明本发明的具体实施。首先，基于文献资料和实测数据进行水动力模型模型初始、边界条件和参数的设定：水动力模型选取渤海海区，采用正交网格，分辨率设置为1/30
°×
1/30
°
。模型的水位、流速开边界条件分别采用nao.99b和tpxo数据，通过对nao.99b和tpxo数据进行克里金插值分别生成模型开边界各个网格点的水位和流速；模型曼宁系数来自现有论文，取为74；模型海底地形来自etopo2数据，波浪模型的风场来自ecmwf数据；采用smagorinsky方程，涡粘系数取为0.28。
19.水动力模型的运算中，将模型设置为冷启动，通过在模型原有曼宁系数的基础上叠加无偏的高斯随机数生成参数集合，使模型自由积分至稳定；采用强化的参数校正数据同化方法（data assimilation scheme for enhancive parameter correction，daepc）进行数据同化，该方法的实现方式为：先进行状态估计，在数值模型达到准平衡状态后，再启动对模型参数的估计。其中，采用eakf方法进行数值模型的状态估计和参数估计。对于本实例，首先基于潮位站测得的潮位、潮流数据和波浪浮漂测得的波高数据对潮位、潮流和波高等状态变量进行数据同化，待模型达到准平衡后，再同时进行状态估计和对曼宁系数的参数估计；上述过程需要反复调试数据同化算法中的局地化半径和膨胀系数（简称同化参数），以获取同化效果最优的模型状态和参数，该过程通过人为设置“真实”解并基于“真实”解生成观测场，测试不同局地化半径和膨胀系数下同化结果与“真实”解的误差，最终确定一组最优的同化参数，本实例中，通过上述理想试验过程对比不同同化参数下模型同化前后的时间、空间平均均方根误差（root mean square error，rmse），最终确定最优局地化半
径为15个网格点，状态膨胀系数为1.03，参数膨胀系数为1.10；基于水动力模型的再分析结果，对水环境模型的水动力部分进行设置，本实例中，水环境模型设置为大连新机场区域（39
°
n-39.4
°
n，121.4
°
e-121.75
°
e）同样采用正交网格，模型分辨率200m
×
200m，曼宁系数采用eakf优化后的结果，采用水动力模型最后一个时刻的曼宁系数集合进行水环境数值模拟，水环境模型的潮汐潮流边界条件采用水动力模型的数值模拟结果。水环境模型温度设置为15
°
c，盐度取为30
‰
，风场仍然采用ecmwf数据。氨氮浓度采用零梯度开边界条件。基于实测氨氮和硝酸盐氮的平均值给定水环境模型初始场。
20.采用与水动力模型相同的方式测试同化参数，最终确定局地化半径为10个网格点，状态膨胀系数为1.06，参数膨胀系数1.05时，模型时、空平均rmse最小，为最优同化参数。
21.使水环境模型自由积分至稳定，同样采用daepc方法，首先基于潮位、潮流、波高和氨氮浓度数据进行水环境模型状态估计，待水环境模型达到准平衡状态，同时进行状态估计和氨氮扩散系数的参数估计。同化过程再次稳定后，保存模型分析结果，最终能够获得高精度、时空连续的氨氮浓度分布数据集。
22.传统的enkf方法直接基于集合成员与集合均值来计算状态变量的分析误差协方差矩阵，以表示同化后状态变量的变化特征：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（1）其中，为分析误差协方差矩阵；为第i个状态向量集合样本；为集合均值；为单位矩阵；为卡尔曼增益矩阵，是集合成员在同化过程中的调整系数；为状态转移矩阵，用于将状态向量投影到观测的状态空间上；为背景误差协方差矩阵，即同化前集合成员之间的协方差。
23.然而，根据kalman滤波的理论推导公式，背景误差协方差理论大小为。显然，公式（1）中的计算方法将使分析误差协方差被低估，传统的enkf通过人为在观测数据上叠加无偏的高斯随机数来提高分析误差协方差的大小，这种方式将导致同化精度降低。而本发明中采用的eakf方法通过设定状态-观测联合向量（）改变了分析误差协方差的计算形式，此时联合向量的背景误差协方差理论公式可变为：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（2）其中，为联合向量的状态转移矩阵，为观测误差协方差矩阵。
24.对上述公式进行奇异值分解，可得到如下形式：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（3）
其中，为一线性算子，可以对联合向量的每个样本进行更新。根据公式（3），能够推导出联合向量集合成员的更新公式：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（4）其中，表示同化前，表示同化后。
25.通过上述过程，eakf消除了传统enkf方法中人为导致的观测误差。此外，由于不需要进行观测的样本化，eakf同化过程的确定性更强，所需的集合样本数量更少，这也说明eakf方法相对传统的enkf需要的计算资源更少，运算速度更快。
26.水动力模型的运算中，将模型设置为冷启动，给定各时刻潮位、潮流数据作为边界条件，曼宁系数设置为一常数，需要说明的是，根据海区具体海底地貌，曼宁系数取值不同，可参照相关文献或已有研究结果设定模型的曼宁系数值。
27.水环境模型的运算中，将模型设置为热启动。
28.实测数据包括地形、风场、氨氮、硝酸盐氮、海区温度、海区盐度等初始场和数据同化观测数据，以及潮位、潮流边界数据。
29.文献资料包括ecwmf的风场资料，以及验潮站的潮位、潮流数据，以及nao.99b和tpxo的潮位、潮流数据。
30.水动力模型的边界条件采用nao.99b潮位数据和tpxo潮流数据。
31.本发明中数据同化模块采用基于最小二乘框架的集合调整卡尔曼滤波方法（eakf），以实测数据对模型的潮位、潮流、波高、曼宁系数、水质浓度和扩散系数进行订正。该方法的实施包括如下两个步骤：第一步，计算观测增量：第一步，计算观测增量：
ꢀꢀꢀꢀꢀ
（5）式中，为第i个集合成员在观测点上的观测增量；为观测误差的标准差；为集合在该观测点上的先验（同化前）误差标准差；为第i个集合成员在该观测点上的先验值；为某一观测数据；为投影到该观测点上的集合平均值。
32.第二步，将观测增量投影到模式网格点上，获得集合成员的更新值：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（6）式中，为第i个先验集合成员；为状态集合与观测之间的协方差；为局地化因子，采用gaspari-cohn函数形式：cohn函数形式：cohn函数形式：
ꢀꢀꢀꢀꢀꢀꢀ
（7）式中，为经验给定的局地化半径；为观测点与状态变量格点之间的距离。
33.为保证集合的离散度，避免滤波发散。对于状态估计，本发明引入静态膨胀方案。即确定一个常数膨胀因子，对各个集合成员相对于集合平均的扰动进行膨胀，用以调整集合的离散度，从而避免集合离散度降低。而对于参数估计，则引入条件静态膨胀方案，即通过判断当前时刻集合方差的大小，决定是否对其进行参数膨胀。条件静态乘法膨胀公式如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（8）式中，为膨胀前、后的参数集合；为参数的集合平均；为根据经验设定的参数膨胀系数；和分别为初始集合和t时刻集合的标准差。
34.基于上述方法，进行数据同化。数据同化过程中，还采用了daepc方法：首先仅进行状态估计至准平衡状态。待模型稳定，加入参数估计并同化至模型再次稳定。
35.本发明创造通过将实测数据与数值模拟结果相结合，提升了水环境数值模拟的精度，能够获得高精度水质再分析数据，为工程决策提供较为准确的再分析数据支撑，有利于进行可行性和安全性的评估；通过改进模型初始场和参数，能够实现短期高精度水环境数值预报，为污水排放、海岸工程等项目提供参考和建议；通过方案设计和集合调整卡尔曼滤波（eakf）方法的使用，提高了参数估计的可靠性和精度，有效减少了计算时间。
36.以上所述仅为本发明创造的较佳实施例而已，并不用以限制本发明创造，凡在本发明创造的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明创造的保护范围之内。