CN109086948A - 基于数据同化的湖库富营养化预警方法 - Google Patents

Abstract

本发明提供一种基于数据同化的湖库富营养化预警方法，包括：获得研究区内的时间上连续的地面监测指标的监测数据；获得研究区的某个时间点监测到的卫星遥感数据；建立研究区的水动力水质水生态过程模型；采用数据同化算法，实现模型结果与监测数据段的同化，水动力水质水生态过程模型进行水质预测，实现湖库富营养化预警。优点为：由于本发明对地面监测数据、卫星遥感数据以及水动力水质水生态过程模型预测到的监测指标预测值这三类数据进行了多源多尺度数据融合，在水动力水质水生态过程模型预测过程中，采用同化算法对水动力水质水生态过程模型存在的不确定性进行修正，从而全面提高了湖库富营养化预警的准确度。

Description

基于数据同化的湖库富营养化预警方法

技术领域

本发明属于湖库富营养化预警技术领域，具体涉及一种基于数据同化的湖库富营养化预警方法。

背景技术

湖库富营养化是指湖库水体接纳过量的氮、磷等营养性物质，使水体中藻类以及其他水生生物异常繁殖，水体透明度和溶解氧变化，造成湖库水体恶化，加速湖库老化，从而使湖库生态和水功能受到阻碍和破坏。对湖库水生态状态进行预警监测，是判断其发展趋势以及制定相应对策的重要手段。

现有的湖库富营养化预警方法，只能依靠单一数据源进行监测预警，从而导致湖库富营养化预警的准确性较低。

发明内容

针对现有技术存在的缺陷，本发明提供一种基于数据同化的湖库富营养化预警方法，可有效解决上述问题。

本发明采用的技术方案如下：

本发明提供一种基于数据同化的湖库富营养化预警方法，包括以下步骤：

步骤1，确定研究区内与湖库富营养化相关的地面监测指标以及确定研究区内的监测点位；

对研究区内的每个监测点位按监测频次进行持续水生态监测，获得时间上连续的地面监测指标的监测数据，所述地面监测指标的监测数据包括气象监测数据、水文监测数据和水质监测数据；假设共有P个时间上连续的地面监测指标的监测数据序列，分别对应监测时刻T₁,T₂,…,T_P；

对研究区进行卫星遥感监测，获得某个时间点监测到的卫星遥感影像；对所述卫星遥感影像进行分析，得到研究区在某时间点所对应的遥感监测指标的时空分布场，遥感监测指标的时空分布场在水体的表层为时间上不连续，空间连续的分布，分别对应监测时刻T_a,T_b,…,T_Q；

步骤2，建立研究区的水动力水质水生态过程模型；其中，将步骤1获得的T₁时刻对应的地面监测指标的监测数据称为初始观测数据；将步骤1获得的T₂,…,T_P时刻对应的地面监测指标的监测数据和T_a,T_b,…,T_Q时刻对应的遥感监测指标的监测数据称为观测数据；

所述水动力水质水生态过程模型用于对步骤1得到的所述地面监测指标的监测数据和所述遥感监测指标的时空分布场进行同化分析，得到研究区的水动力、水质和水生态的变化规律；

然后，所述水动力水质水生态过程模型采用以下方式运行：

步骤2.1，预测过程：所述水动力水质水生态过程模型以T₁时刻对应的所述初始观测数据作为初始场，使用在计算区域边界的地面监测指标的监测数据驱动所述水动力水质水生态过程模型运算，预测得到之后任意某个时刻所对应的监测指标预测值；

步骤2.2，如果步骤2.1得到的任意某个时刻所对应的监测指标预测值为T₂时刻的监测指标预测值，则对T₂时刻的监测指标预测值和T₂时刻所对应的地面监测指标的监测数据进行加权分析，得到T₂时刻最优模型驱动的边界输入值，并采用所述边界输入值驱动所述水动力水质水生态过程模型运算，预测得到之后任意某个时刻所对应的监测指标预测值；

步骤2.3，如果步骤2.2预测得到的任意某个时刻所对应的监测指标预测值为T_a时刻的监测指标预测值；

则采用T_a时刻的遥感监测指标的监测数据，对T_a时刻的监测指标预测值进行同化计算，得到新的T_a时刻的监测指标预测值；然后，采用所述新的T_a时刻的监测指标预测值替代原T_a时刻的监测指标预测值，成为步骤2.1中的新初始场，不断向前积分直到预测得到T_b时刻所对应的监测指标预测值；

如此不断迭代计算，直到完成所有有观测数据时刻的状态预测和更新，得到最终优化后的水动力水质水生态过程模型；

步骤3，采用步骤2最终得到的水动力水质水生态过程模型对未来某个时间的监测指标时空分布场进行预测，并判断预测到的监测指标时空分布场是否在设定阈值内；如果在，则表明监测指标正常，不进行预警；如果不在，则表明监测指标异常，进行湖库富营养化预警。

优选的，所述监测指标包括水动力监测指标、水温监测指标和水质监测指标。

优选的，所述水质监测指标包括叶绿素、总磷、总氮、高锰酸盐指数和溶解氧。

优选的，步骤2具体为：

步骤2-1，所述水动力水质水生态过程模型以T₁时刻对应的所述初始观测数据作为初始场，不断向前积分计算，得到在某个监测点和某个监测时间对应的监测指标预测值；

步骤2-2，判断是否存在与某个监测点和某个监测时间对应的遥感监测指标观测值；如果存在，则执行步骤2-3至步骤2-7；如果不存在，则执行步骤2-6至步骤2-7；

步骤2-3，按下式计算卡尔曼增益矩阵K：

其中：

s为单元号；l为水体某层；

K(s,l)：卡尔曼增益矩阵；

P_PL(s,l)：下层预测变量的方差；

C(s,l)：预测值到观测值之间的线性映射的参数；

C^T(s,l)：预测值到观测值之间的线性映射的参数的转置矩阵

R(s,l)：：与多尺度相关的观测值的噪声方差；

步骤2-4，按下式将监测指标预测值映射到观测空间：

将预测值映射到观测空间的观测值的估计值；

由下预测到本层状态变量x的均值；

D(s,l)：预测值到观测值之间的线性映射的参数；

新息表达式为：

IN(s,l)：新息；

y(s,l)：观测值；

观测值的估计值；

按下式输出监测指标更新估计值：

监测指标的更新估计值；

由下预测到本层状态变量x的均值；

K(s,l)：卡尔曼增益矩阵；

步骤2-5，按下式计算监测指标更新估计值的方差；然后执行步骤2-7；

P_U(s,l)＝[I–K(s,l)C(s,l)]P_PL(s,l)

P_U(s,l)：更新估计值的方差；

I为m×m维的恒等矩阵，即单元矩阵；

步骤2-6，监测指标更新估计值直接表示为监测指标预测值，然后执行步骤2-7；

步骤2-7，将更新得到的监测指标更新估计值作为初始场，更新步骤2-1的初始场，对水动力水质水生态过程模型的初始参数进行更新，由此迭代循环。

本发明提供的基于数据同化的湖库富营养化预警方法具有以下优点：

由于本发明对地面监测数据、卫星遥感数据以及水动力水质水生态过程模型预测到的监测指标预测值这三类数据进行了多源数据融合，并在水动力水质水生态过程模型每一步的预测过程中，如果存在卫星遥感数据，则采用卫星遥感数据及时对水动力水质水生态过程模型偏差进行修正，从而全面提高了湖库富营养化预警的准确度。

附图说明

图1为本发明提供的基于数据同化的湖库富营养化预警系统的功能结构图；

图2为本发明提供的数据同化过程的流程图；

图3为本发明提供的空间滤波和平滑过程的流程图。

具体实施方式

为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

湖库通常为重要的饮用水水源地，湖库富营养化问题更是严重威胁着人类的饮用水安全。湖库水华及富营养化的发生受气象、水动力学条件、污染负荷输入、湖库水质等的共同影响。

本发明中，进行库富营养化预警的数据来源共有三种：1.卫星遥感数据；2.地表监测数据；3.水动力水质水生态过程模型预测到的模型数据，该模型数据即是水动力水质水生态过程模型预测到的监测指标预测值。

具体的，水体因为各组分及其含量不同造成水体吸收和散射的变化，使一定波长范围反射率显著不同，是定量估测内陆水体水质参数的基础，也会是湖泊富营养化遥感的理论基础。水体中影响光谱反射率的物质主要有三类：(1)浮游植物，主要是各种藻类；(2)由浮游物质死亡而产生的有机碎屑及泥沙等无极悬浮颗粒，总称悬浮物；(3)由黄腐酸、腐殖酸等组成的溶解性有机物。除了悬浮物在自然浓度条件下对光不发生明显的吸收外，其余两种物质分别选择吸收一定波长范围的光，形成各自的特征吸收波谱。同时这些物质对光的散射使得光改变方向，其中后项散射光与水底的反射光一起返回水面，通过水气界面回到大气中，是可以遥测到的部分。光活性物质包括浮游藻类色素(叶绿素)、非色素悬浮颗粒以及有色可溶性有机物(CDOM)，藻类色素是决定水体反射光谱特征的重要因素，也是水色遥感的重要参数。

湖库富营养化特征之一是藻类生长和聚集。利用多光谱数据遥感，对叶绿素进行反演，得到湖库在某一时刻上湖库表层空间上的叶绿素分布，得到富营养化状态。卫星遥感数据反演了湖库的上层的叶绿素分布水平和状态，且仅为叶绿素等有限指标。卫星遥感数据，受卫星过境的周期和当时的天气如云等的限制，呈现时间上间断，空间上较为连续的数据结果。

地面时间监测站，受区域限制，通常在特定点位有间隔的数据输出，通常气象和水文数据间隔为1天，水质监测的指标通常根据《地表水环境质量标准》(GB3838-2002)中24项指标进行监测，通常在水文站或关键控制站点布设站点，一般水质站点丰、平、枯水期各两次，重要站点每月1次。受时间和空间点位限制，地面观测数据通常不能捕捉到局部水体的水华暴发，但是地面观测数据能得到水体较为全面的状态指标，也是驱动水动力水质水生态过程模型的基础及关键变量。

为准确预报湖库的富养养化状态，仅用的卫星遥感及地面监测数据都难以实现。本发明主要思路为：本发明以地表监测数据驱动水动力水质水生态过程模型运算，结合任意M1时刻的气象监测数据、水文监测数据和水质监测数据等，得到任意M1时刻的监测指标预测值；然后，如果M1时刻存在卫星遥感数据，则采用M1时刻的卫星遥感数据对M1时刻的监测指标预测值进行同化计算，得到一个新的M1时刻的监测指标预测值；

然后，以新的M1时刻的监测指标预测值作为初始场，驱动水动力水质水生态过程模型运算，在水动力水质水生态过程模型运算过程中，需要结合之后任意M2时刻的气象监测数据、水文监测数据和水质监测数据等，从而预测得到M2时刻的监测指标预测值；如此不断向前迭代运算，实现对准确预报湖库的富养养化状态。

因此，本发明采用同化算法对卫星遥感数据、地表监测数据、和水动力水质水生态过程模型预测到的模型数据，进行多源数据融合，同时水动力水质水生态过程模型给出空间和时间连续的监测指标值，对湖库富营养化进行预警。

具体的，参考图1和图2，本发明提供的一种基于数据同化的湖库富营养化预警方法，包括以下步骤：

步骤1，确定研究区内与湖库富营养化相关的地面监测指标以及确定研究区内的监测点位；其中，监测指标包括但不限于水动力监测指标、水温监测指标和水质监测指标。水质监测指标包括但不限于叶绿素、总磷、总氮、高锰酸盐指数和溶解氧。

对研究区内的每个监测点位按监测频次进行持续水生态监测，获得时间上连续的地面监测指标的监测数据，所述地面监测指标的监测数据包括气象监测数据、水文监测数据和水质监测数据；假设共有P个时间上连续的地面监测指标的监测数据序列，分别对应监测时刻T₁,T₂,…,T_P；

对研究区进行卫星遥感监测，获得某个时间点监测到的卫星遥感影像；对所述卫星遥感影像进行分析，得到研究区在某时间点所对应的遥感监测指标的时空分布场，遥感监测指标的时空分布场在水体的表层为时间上不连续，空间连续的分布，分别对应监测时刻T_a,T_b,…,T_Q；

步骤2，建立研究区的水动力水质水生态过程模型；其中，将步骤1获得的T₁时刻对应的地面监测指标的监测数据称为初始观测数据；将步骤1获得的T₂,…,T_P时刻对应的地面监测指标的监测数据和T_a,T_b,…,T_Q时刻对应的遥感监测指标的监测数据称为观测数据；

所述水动力水质水生态过程模型用于对步骤1得到的所述地面监测指标的监测数据和所述遥感监测指标的时空分布场进行同化分析，得到研究区的水动力、水质和水生态的变化规律；

实际应用中，水动力水质水生态过程模型是一个剖面二维(纵向和垂向)、横向平均的水动力学水质水生态过程模型，可用来模拟河流、河口、水库和流域系统。由于该模型呈河宽方向均匀，因此最适用于模拟展示纵向和垂向水质梯度的相对狭长的水体，该模型同时也适合模拟一些具有湖泊特性的河流。模型的富营养化过程涉及温度－营养物－藻类－溶解氧－有机物－沉积物之间的关系，可模拟包括DO、BOD、大肠杆菌、藻类等在内的21种水质变量浓度变化。模型可以模拟有坡度的、受潮汐影响的、密度分层的河流的二维水动力和水质，且有最新的模型改进，减少了数值计算的错误，提高了模拟精度。

水动力水质水生态过程模型基本方程包括：

连续方程

X方向的动量方程

Z方向的动量方程

水位波动方程

紊流方程

式中：为紊动动能生成项；

为浮力生成项，该浮力项在稳定分层时可抑制紊动动能的生成，削弱热量向下的传递，是水库能保持稳定分层的重要因素；

ν_e——分子粘性系数ν与紊动涡粘系数ν_t之和，

m2/s，U、W——纵向和垂向流速，

m/s；p——压强，Pa；

B——水体宽度，为高程的函数，m；

K——紊动动能，m²/s²；

Ε-紊动动能耗散率，m2/s3；

σ_k、σ_ε——分别为紊动动能和耗散率的普朗特数，一般取1.0和1.3。

C_μ、C_ε1、C_ε2——模型常数，取值分别为0.09、1.44、1.92。

水动力学模型方程组中的常数为通用常数，其取值由基本实验确定，不因具体问题而改变。

C_ε3——为铅垂方向速度w和水平方向速度u比值的函数，取值介于0～1。水温基本方程为：

式中：

T-水温，℃；

B-宽度，m；

u-水平方向的流速分量，m/s；

w-垂直方向的流速分量，m/s；

D_tx-横向的水温紊动扩散系数；

D_tz-垂向的水温紊动扩散系数；

φ-太阳热辐射通量，J/(m²·s)；

ρ-水体密度，kg/m³；

C_p-水的比热，J/(kg·℃)；S(T)-温度源项；

同理，水质数学模型的基本方程为：

式中：

C-水质指标浓度，mg/L；

S(C)—污染物源(汇)项；

f(C)—污染物反应项；

本水动力水质水生态过程模型可以模拟以下参数：水位、流速、水温、溶解氧、大肠杆菌、无极悬浮物、浮游植物、CBOD、氨氮、硝氮-亚硝氮、磷酸盐、溶解态有机物、颗粒态有机物、总无机碳、碱度、总铁、有机沉积物。

本水动力水质水生态过程模型具有以下功能：

(1)水动力学：该模型能够预测水面高程、速度和温度。由于温度对水中各物质浓度有影响,它包含在水动力学计算中。

(2)水质：在水质模拟中可以加入或排除任何组合中的一项或多项成分。由于盐度或全部溶解固体(盐度)对密度产生的影响，只有它们需要被模拟时才加入到水动力学计算中。水质算法模块允许成分能够根据需求增加一项，作为附加的子程序。水质状态变量除了温度之外还包含：保守示踪剂、水的停留时间、大肠杆菌、污染物、无机悬浮固体、浮游植物、植物病原体、CBOD、铵、硝酸盐和亚硝酸盐、生物利用的磷(一般由正磷酸盐和可溶的磷代表)、易溶解的有机物、难溶的有机物、易降解的颗粒有机物、难降解的颗粒有机物、总无机碳、碱浓度、总铁、溶解氧、有机沉积物、气体挟带等。此外，模型可计算出60多个派生变量包括pH值、TOC、DOC、TON、TOP、DOP等等内部的状态变量并输出，以与实测数据对比。

(3)多个分支：分支算法可以应用在形状复杂的水体，如树枝状的水库或河口。

(4)多个水体：该模型可用于任何数量的河流、湖泊、水库、河口及其组合。

(5)可变的网格间距：变量段的长度和层高允许需要的高分辨率的规格，水体之间可以采用不同的垂向网格间距。

(6)水质不依赖于水动力学：水质可以相对水动力更低的频率更新，从而减少对计算的要求。然而，水质与水动力是分离的(即分开、独立的水动力和水质代码，从水动力模型中的输出是储存在硬盘上，然后指定用于水质计算的平流流量)。对于长期的水动力输出用来驱动水质模型的存储要求被禁止，除了非常小的网格。此外，当用来驱动水质的水动力数据输入每个步长时，减少的计算时间是最小的。

(7)自适应步长：该模型包含一个变量的步长算法，以确保水动力学所强加的数值解方案不受破坏，从而达到稳定要求。

然后，所述水动力水质水生态过程模型采用以下方式运行：

步骤2.1，预测过程：所述水动力水质水生态过程模型以T₁时刻对应的所述初始观测数据作为初始场，使用在计算区域边界的地面监测指标的监测数据驱动所述水动力水质水生态过程模型运算，预测得到之后任意某个时刻所对应的监测指标预测值；

步骤2.2，如果步骤2.1得到的任意某个时刻所对应的监测指标预测值为T₂时刻的监测指标预测值，则对T₂时刻的监测指标预测值和T₂时刻所对应的地面监测指标的监测数据进行加权分析，得到T₂时刻最优模型驱动的边界输入值，并采用所述边界输入值驱动所述水动力水质水生态过程模型运算，预测得到之后任意某个时刻所对应的监测指标预测值；

步骤2.3，如果步骤2.2预测得到的任意某个时刻所对应的监测指标预测值为T_a时刻的监测指标预测值；

则采用T_a时刻的遥感监测指标的监测数据，对T_a时刻的监测指标预测值进行同化计算，得到新的T_a时刻的监测指标预测值；然后，采用所述新的T_a时刻的监测指标预测值替代原T_a时刻的监测指标预测值，成为步骤2.1中的新初始场，不断向前积分直到预测得到T_b时刻所对应的监测指标预测值；

如此不断迭代计算，直到完成所有有观测数据时刻的状态预测和更新，得到最终优化后的水动力水质水生态过程模型；

步骤2具体为：

步骤2-1，所述水动力水质水生态过程模型以T₁时刻对应的所述初始观测数据作为初始场，不断向前积分计算，得到在某个监测点和某个监测时间对应的监测指标预测值；

步骤2-2，判断是否存在与某个监测点和某个监测时间对应的遥感监测指标观测值；如果存在，则执行步骤2-3至步骤2-7；如果不存在，则执行步骤2-6至步骤2-7；

步骤2-3，按下式计算卡尔曼增益矩阵K：

其中：

s为单元号；l为水体某层；

K(s,l)：卡尔曼增益矩阵；

P_PL(s,l)：下层预测变量的方差；

C(s,l)：预测值到观测值之间的线性映射的参数；

C^T(s,l)：预测值到观测值之间的线性映射的参数的转置矩阵

R(s,l)：：与多尺度相关的观测值的噪声方差；

卡尔曼增益矩阵K表示观测值对预测值的影响的程度。观测值方差越高(与预测值方差相比)，影响越小。卡尔曼增益矩阵K的每个值的范围都在0到1之间，0表示观测值的影响可以被忽略；1表示观测值为唯一重要的估计量。

步骤2-4，按下式将监测指标预测值映射到观测空间：

将预测值映射到观测空间的观测值的估计值；

由下预测到本层状态变量x的均值；

D(s,l)：预测值到观测值之间的线性映射的参数；

新息表达式为：

IN(s,l)：新息；

y(s,l)：观测值；

观测值的估计值；

新息(innovation)，表示观测值的估计值与实际观测值之间的差异，也即更新值随观测值的变化程度。

按下式输出监测指标更新估计值：

监测指标的更新估计值；

由下预测到本层状态变量x的均值；

K(s,l)：卡尔曼增益矩阵；

步骤2-5，按下式计算监测指标更新估计值的方差；然后执行步骤2-7；

P_U(s,l)＝[I–K(s,l)C(s,l)]P_PL(s,l)

P_U(s,l)：更新估计值的方差；

I为m×m维的恒等矩阵，即单元矩阵；

步骤2-6，监测指标更新估计值直接表示为监测指标预测值，然后执行步骤2-7；

式中，结合了当前层观测信息和下层变量的估计。

步骤2-7，将更新得到的监测指标更新估计值作为初始场，更新步骤2-1的初始场，对水动力水质水生态过程模型的初始参数进行更新，由此迭代循环。

步骤3，采用步骤2最终得到的水动力水质水生态过程模型对未来某个时间的监测指标时空分布场进行预测，并判断预测到的监测指标时空分布场是否在设定阈值内；如果在，则表明监测指标正常，不进行预警；如果不在，则表明监测指标异常，进行湖库富营养化预警。

本发明提供的基于数据同化的湖库富营养化预警方法，在具体实现上，通过基于数据同化的湖库富营养化预警平台实现，湖库富营养化预警平台分数据层、数据同化和水动力水质水生态过程模型三大功能部分，其中水动力水质水生态过程模型包含了过程模型和举动模型的数据输入和反映变量。

实际应用中，可以定义数据同化的参量为叶绿素浓度，利用数据同化算法实现水动力水质水生态过程模型、地面监测数据和卫星遥感数据的同化处理。湖库富营养化预警平台运行后，用户通过交互式对话定义输入和输出文件的格式、路径和文件，确定同化叶绿素时间空间分布，并显示流程处理状态。

在仪器观测和模型模拟分别有不间的空间尺度时，应用多尺度卡尔曼滤波(MKS-EM)框架或集合卡尔曼滤波算法(EnKF)，可以很好地提高状态变量的预测精度，同时框架还可以用于数据融合。MKS框架中需要的参数包括① 不同尺度上的所有观测或模型预测值的平均值；②P(root,1)；③Q(m,l)；④R(m,l)。其中，后三者都是必要的输入。

基于MKS框架可总结为以下步骤：

1.输入参数的初始化。包括4个参数：湖库顶层最大尺度叶绿素估计值的无条件平均值，可看作是所有不同尺度模型预测值的平均值或是直接输入值；P(root,1)相遇无条件预测方差；Q(m,l)多尺度或大尺度动态方程噪声方差；R(m,l)：模型方程噪声方差。

2.数据更新

对于全部水体的变量多尺度分布，首先采用输入的条件，对水体进行以下变量的预测，包括不同深度及时间上的水动力学、水温、水质、叶绿素等，如果某一时刻存在卫星遥感影像分析出的水体表层叶绿素观测值，结合模型的预测，进行卡尔曼滤波，更新叶绿素状态变量，如果没有卫星遥感影像分析出的水体表层叶绿素观测值，则将模型的本次预测值直接作为下一步的更新值。

在同化开始时，必须根据观测信息给定这些参数的初始值。基于参数的初始值，MKS算法开始同化计算。同化完毕后，则进入EM算法的参数更新步骤，即计算完全统计量，更新初始参数。在此过程中，需要计算隐藏状态变量的似然函数值。在更新完毕后，则进入下一轮的MKS算法同化计算，如此重复各步骤，直至似然函数值不变化或者变化范围小于某一个设定值。此时，估计得到的隐藏状态变量的数值被认为是最优值的。

驱动湖库水体的富营养化预警预报模型，不仅需要水文、气象等的过程的输入，还需要水体的水质及营养状态过程，通过地面监测数据驱动水动力水质水生态过程模型运行，同时采用卫星遥感数据对水动力水质水生态过程模型的预测结果进行比对，采用数据同化算法(卡尔曼滤波等)，实现数据融合，提高水动力水质水生态过程模型预测精度。

传统的湖库富营养化预警方法，主要存在以下不足：(1)数据来源单一，无法对不同来源的数据进行融合预报，导致富营养化预警的准确度低；(2)无法对水动力水质水生态过程模型的不确定性进行及时纠正，导致在迭代过程中，模型偏差会不断传递，从而降低了富营养化预警的准确度。

而本发明提供的基于数据同化的湖库富营养化预警方法，采用气象、水文、污染负荷等地面监测数据驱动水动力水质水生态过程模型的运算，产生空间上和时间上连续的监测指标预测值，包括水动力、水温、10-20个水质指标的时空分布；然后，采用数据同化算法，将水动力水质水生态过程模型预测到的监测指标预测值，与同一时刻的卫星遥感数据进行数据同化，使得在某个空间和时间上水动力水质水生态过程模型的预测结果逼近实际观测数据，修正初始场，继续水动力水质水生态过程模型的下一步预测，达到对连续时间和空间上湖库富营养化的预警预测。可见，由于本发明对地面监测数据、卫星遥感数据以及水动力水质水生态过程模型预测到的监测指标预测值这三类数据进行了多源多尺度数据融合，并在水动力水质水生态过程模型每一步的预测过程中，如果存在卫星遥感数据，则采用卫星遥感数据及时对水动力水质水生态过程模型偏差进行修正，从而全面提高了湖库富营养化预警的准确度。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (4)

1.一种基于数据同化的湖库富营养化预警方法，其特征在于，包括以下步骤：

步骤1，确定研究区内与湖库富营养化相关的地面监测指标以及确定研究区内的监测点位；

对研究区内的每个监测点位按监测频次进行持续水生态监测，获得时间上连续的地面监测指标的监测数据，所述地面监测指标的监测数据包括气象监测数据、水文监测数据和水质监测数据；假设共有P个时间上连续的地面监测指标的监测数据序列，分别对应监测时刻T₁,T₂,…,T_P；

对研究区进行卫星遥感监测，获得某个时间点监测到的卫星遥感影像；对所述卫星遥感影像进行分析，得到研究区在某时间点所对应的遥感监测指标的时空分布场，遥感监测指标的时空分布场在水体的表层为时间上不连续，空间连续的分布，分别对应监测时刻T_a,T_b,…,T_Q；

步骤2，建立研究区的水动力水质水生态过程模型；其中，将步骤1获得的T₁时刻对应的地面监测指标的监测数据称为初始观测数据；将步骤1获得的T₂,…,T_P时刻对应的地面监测指标的监测数据和T_a,T_b,…,T_Q时刻对应的遥感监测指标的监测数据称为观测数据；

所述水动力水质水生态过程模型用于对步骤1得到的所述地面监测指标的监测数据和所述遥感监测指标的时空分布场进行同化分析，得到研究区的水动力、水质和水生态的变化规律；

然后，所述水动力水质水生态过程模型采用以下方式运行：

步骤2.1，预测过程：所述水动力水质水生态过程模型以T₁时刻对应的所述初始观测数据作为初始场，使用在计算区域边界的地面监测指标的监测数据驱动所述水动力水质水生态过程模型运算，预测得到之后任意某个时刻所对应的监测指标预测值；

步骤2.2，如果步骤2.1得到的任意某个时刻所对应的监测指标预测值为T₂时刻的监测指标预测值，则对T₂时刻的监测指标预测值和T₂时刻所对应的地面监测指标的监测数据进行加权分析，得到T₂时刻最优模型驱动的边界输入值，并采用所述边界输入值驱动所述水动力水质水生态过程模型运算，预测得到之后任意某个时刻所对应的监测指标预测值；

步骤2.3，如果步骤2.2预测得到的任意某个时刻所对应的监测指标预测值为T_a时刻的监测指标预测值；

则采用T_a时刻的遥感监测指标的监测数据，对T_a时刻的监测指标预测值进行同化计算，得到新的T_a时刻的监测指标预测值；然后，采用所述新的T_a时刻的监测指标预测值替代原T_a时刻的监测指标预测值，成为步骤2.1中的新初始场，不断向前积分直到预测得到T_b时刻所对应的监测指标预测值；

如此不断迭代计算，直到完成所有有观测数据时刻的状态预测和更新，得到最终优化后的水动力水质水生态过程模型；

步骤3，采用步骤2最终得到的水动力水质水生态过程模型对未来某个时间的监测指标时空分布场进行预测，并判断预测到的监测指标时空分布场是否在设定阈值内；如果在，则表明监测指标正常，不进行预警；如果不在，则表明监测指标异常，进行湖库富营养化预警。

2.根据权利要求1所述的基于数据同化的湖库富营养化预警方法，其特征在于，所述监测指标包括水动力监测指标、水温监测指标和水质监测指标。

3.根据权利要求2所述的基于数据同化的湖库富营养化预警方法，其特征在于，所述水质监测指标包括叶绿素、总磷、总氮、高锰酸盐指数和溶解氧。

4.根据权利要求1所述的基于数据同化的湖库富营养化预警方法，其特征在于，步骤2具体为：

步骤2-1，所述水动力水质水生态过程模型以T₁时刻对应的所述初始观测数据作为初始场，不断向前积分计算，得到在某个监测点和某个监测时间对应的监测指标预测值；

步骤2-2，判断是否存在与某个监测点和某个监测时间对应的遥感监测指标观测值；如果存在，则执行步骤2-3至步骤2-7；如果不存在，则执行步骤2-6至步骤2-7；

步骤2-3，按下式计算卡尔曼增益矩阵K：

其中：

s为单元号；l为水体某层；

K(s,l)：卡尔曼增益矩阵；

P_PL(s,l)：下层预测变量的方差；

C(s,l)：预测值到观测值之间的线性映射的参数；

C^T(s,l)：预测值到观测值之间的线性映射的参数的转置矩阵R(s,l)：：与多尺度相关的观测值的噪声方差；

步骤2-4，按下式将监测指标预测值映射到观测空间：

将预测值映射到观测空间的观测值的估计值；

由下预测到本层状态变量x的均值；

D(s,l)：预测值到观测值之间的线性映射的参数；

新息表达式为：

IN(s,l)：新息；

y(s,l)：观测值；

观测值的估计值；

按下式输出监测指标更新估计值：

监测指标的更新估计值；

由下预测到本层状态变量x的均值；

K(s,l)：卡尔曼增益矩阵；

步骤2-5，按下式计算监测指标更新估计值的方差；然后执行步骤2-7；

P_U(s,l)＝[I–K(s,l)C(s,l)]P_PL(s,l)

P_U(s,l)：更新估计值的方差；

I为m×m维的恒等矩阵，即单元矩阵；

步骤2-6，监测指标更新估计值直接表示为监测指标预测值，然后执行步骤2-7；

步骤2-7，将更新得到的监测指标更新估计值作为初始场，更新步骤2-1的初始场，对水动力水质水生态过程模型的初始参数进行更新，由此迭代循环。