湖泊藻类动态模型数据同化模式的改进

2021-08-14何欣霞何梦男陈求稳

水资源保护 2021年4期

李港，陈诚，何欣霞,2，何梦男，邓跃，陈求稳,2

(1.南京水利科学研究院生态环境研究所，江苏南京 210029； 2.重庆交通大学河海学院，重庆 400074；3.中建四局第五建筑工程有限公司，广东深圳 518000)

湖泊是重要的淡水资源库，近些年来，由于营养盐的大量排放以及气候变化的影响，湖泊富营养化及水华问题日益严重[1]。据2018年中国生态环境状况公报显示，我国约有1/3的湖库达到了富营养化状态[2]，其中太湖、滇池、巢湖等重要水源地富营养化形势严峻，严重威胁饮用水安全和人体健康[3]。水华的精确模拟预测是富营养化湖泊风险管控的基础，对于富营养化湖泊的治理具有重要意义。

富营养化模型是进行湖泊水华模拟预测的重要工具。国内外学者开发了EFDC(environmental fluid dynamics code)、WASP(water quality analysis simulation program)、Delft3D-Eco和三维水生态动力学(three-dimensional hydro-ecological dynamics, 3DHED)模型等一系列湖泊水质模型[4-7]，已成功应用于滇池[8]、Langelmavesi-Roine湖[9]、维多利亚湖[10]等富营养化湖泊的藻类模拟预测。但由于大范围蓝藻生物量的空间分布获取困难，模型使用站点观测数据作为输入数据时，会产生较大的初值误差，同时，模型运行过程中模拟误差也会不断积累，使得模拟精度无法得到有效保证[11]。数据同化(data assimilation, DA)技术通过将观测数据融入模型计算过程，校正模型模拟结果，同步更新模型参数，降低了输入数据和模型结构不确定性的影响[12]。该技术包含连续数据同化和顺序数据同化两大类，其中连续顺序同化算法利用同化时间窗口内的所有观测数据和模型状态值进行最优估计，不断调整模型初始场，最终将模型轨迹拟合至同化窗口周期内的观测值上，主要以变分算法为主[13]。相较于变分算法，以集合卡尔曼滤波为代表的顺序数据同化算法能够提供状态量的均值及其相应的误差协方差，其预报误差也可以随着模型动态变化，同时由于引入了集合的思想，不需要伴随或线性算子，具有适用于非线性系统、程序设计相对简单以及易实现并行运算等特点[14]，近年来得到了广泛的应用。Loos等[15]结合HSFP、EFDC以及WQFS(water quality forecast system)，建立了基于EnKF (ensemble Kalman filter)的Yeongsan河藻类数值模型和数据同化系统。刘卓等[16]采用2009—2011年太湖站点监测数据对Delft3D-BLOOM模型进行EnKF数据同化，提高了Delft3D-BLOOM模型对太湖蓝藻的模拟精度。对于大型湖泊而言，数据同化技术虽然能将初始场的使用时间减少到一个观测周期内，降低了模型误差的累积，但加入的同化数据是稀疏站点数据插值得到的蓝藻生物量数据，忽略了蓝藻生物量的空间异质性，存在较大的初值误差。遥感数据由于具有较好的空间代表性[17]，尤其对大型湖区的蓝藻生物量表征能力优于站点观测数据，所以将其与数据同化技术相结合，以降低初值误差，有效提高模型模拟与预测精度。但某些时刻遥感数据质量受云雾等气象条件的影响，观测误差较大，此时加入低质量的遥感数据进行数据同化，反而会造成一个观测周期内的模拟精度降低。

基于此，本文以太湖为研究区域，采用2014—2016年水环境监测数据率定的3DHED模型进行藻类时空模拟，加入蓝藻生物量遥感反演数据进行模型的EnKF数据同化，并提出一种数据同化改进策略，通过比较两种方式模拟结果与实测数据的均方根误差(RMSE)，进行数据同化时遥感数据的有效筛选，避免因部分时刻遥感数据误差较大给模型带来的影响，提升蓝藻生物量模拟精度，以期为富营养化湖泊水华的精确模拟预测和风险管控提供技术支撑。

1 研究区概况与数据资料

1.1 研究区概况

太湖是中国第三大淡水湖泊，总面积为 2 338 km2，平均水深约2 m[18]，分为梅梁湾、竺山湖、贡湖、西北湖、湖心区、东部湖区、西南湖和东太湖8个湖区(图1)。1980年以来，随着太湖周边城市化进程的不断加快，湖泊的不合理利用开发加剧[19]，如工业生产、围湖造田以及沿岸渔业发展等，导致大量污染物汇入太湖，使得太湖生态系统逐步退化，水体自净能力下降，水质日趋恶化，蓝藻水华频发[20]。2007年蓝藻水华暴发造成了无锡市近百万人的饮用水危机，引起全社会广泛关注[21]。2017年，太湖再次暴发大面积蓝藻水华，严重影响周边群众生产生活[22]，太湖富营养化形势依然严峻。

图1 研究区域及采样点分布

1.2 数据资料

2 模型与研究方法

2.1 3DHED模型

3DHED模型通过耦合SELFE模型与SALMO模型，考虑了动力条件对藻类分布的影响，可以有效模拟藻类时空变化[7]。其耦合原理为：①通过SELFE模型将湖泊划分为若干个三角网格，并将SALMO模型应用到其中的每个网格；②通过SELFE模型计算得到全湖所有网格的水动力数据，然后将其提供给SALMO模型；③根据水动力数据与模型水质初始场信息，利用SALMO模型在每个网格上进行水质模拟。3DHED模型参数率定参考郭静等[24]的前期研究成果，采用2014—2016年的水环境生态监测数据对模型影响较大的20个参数进行率定。选定的参数及其率定结果如表1所示。

2.2 EnKF数据同化

EnKF计算过程包括预测和分析两个步骤：设有N个集合，在k=0时刻对每个集合进行初始化。式(1)为预测阶段，式(2)～(6)为分析阶段。

(1)

式中：Xai,k为k时刻第i集合状态变量分析值；Xfi,k+1为k+1时刻状态变量预测值；Mk,k+1为模型算子，本文选取3DHED模型；wi,k为模型结构不确定性引起的模拟误差，服从均值为0、协方差矩阵为Qk的正态分布。

Xai,k+1=Xfi,k+1+Kk+1[Y0k+1-Hk+1(Xfk+1)]

(2)

(3)

表1 模型参数率定结果

(4)

(5)

(6)

同时进行状态变量和模型参数同化，对模型模拟精度的提高有明显效果[25]。本文以太湖优势藻种蓝藻生物量作为状态变量，以敏感性最高的蓝藻最佳生长温度作为模型参数[16]，同时进行状态变量和模型参数数据同化。采用EnKF算法进行数据同化时，合理的同化方案和合适的参数设置是保证同化效果的前提。EnKF算法利用集合思想解决了实际应用中背景协方差矩阵估计和预报困难的问题，集合数越多，集合均值协方差与真值协方差越接近，但计算负荷会显著增加；集合数过少，则可能导致模型误差协方差估计错误。本研究将集合数设置为30、50、70、100、200和500进行数值试验，结果发现集合数为100时可较好地兼顾模拟精度和计算效率。采用数据同化方法将观测数据与模拟数据相融合，发现融合的优劣性与模拟误差协方差和观测误差协方差紧密相关[26]，所以选择合适的模拟误差和观测误差对提高同化精度至关重要。本文参考Chen等[27-28]的研究结果，将观测误差与模拟误差分别设置为1%、10%、20%、30%，筛选观测误差和模拟误差的最优组合，进行同化参数设置。

2.3 改进数据同化

同化模型中背景场和观测数据的误差均会导致模型的预报误差，而背景场与观测误差均与实测数据相关，因此本研究提出一种改进的数据同化(modified data assimilation, mDA)策略。首先在有遥感数据的时刻，加入遥感数据进行模型的数据同化计算，同时进行模型模拟计算，直到有站点实测数据的下一时刻，以站点实测数据为参考，计算并比较该时刻同化结果和模型模拟结果的RMSE，从而进行同化时刻数据源的筛选。提出的mDA能够有效剔除低质量的遥感数据，避免遥感数据质量较差所导致的同化结果精度降低的问题，从而提高同化模型的鲁棒性。改进数据同化流程如图2所示。

图2 改进的数据同化流程

2.4 模型精度评价指标

本文采用RMSE和一致性系数IOA(index of agreement)定量评价蓝藻生物量的观测数据与模拟结果的拟合度。RMSE表示模拟结果与站点实测数据的偏差，其值越小说明模拟结果与站点实测数据越接近。IOA表示模拟结果与站点实测站点数据的吻合程度，描述了模拟结果与站点实测数据在变化趋势上的一致性，IOA值介于0～1之间，越接近1表示模拟结果与实测数据吻合度越高。

(7)

(8)

3 结果与分析

3.1 数据同化参数确定

表2为16种不同观测误差和模拟误差组合下模拟结果与观测数据的RMSE值。从表2中可以看出，当观测误差设置为1%且模拟误差设置为10%时，数据同化结果的RMSE最小，因此本文以观测误差1%、模拟误差10%作为EnKF数据同化方法的参数设置。

表2 不同观测误差与模拟误差下的RMSE值

3.2 数据同化结果分析

针对太湖8个湖区各湖区选择一个代表性站点进行3DHED模拟结果和DA模拟结果的比较，图3为8个代表性站点蓝藻生物量观测值、3DHED模拟结果以及DA模拟结果的对比。从图3可以看出，3DHED模拟结果虽能在一定程度上反映蓝藻生物量的变化趋势，但对蓝藻生物量峰值的捕捉能力欠佳，而加入遥感数据进行同化后，同化时刻初值误差有所降低，模型误差累积周期缩短，使得DA模拟结果的蓝藻生物量变化趋势与峰值的模拟效果均较3DHED模拟结果有显著提高。8个代表性站点的3DHED模拟结果、DA模拟结果与站点实测数据的RMSE和IOA及其平均值见表3。从表3可以看出，8个代表性站点3DHED模拟结果的RMSE均值为 1.44 mg/L，IOA均值为0.43；DA模拟结果的RMSE均值为1.34 mg/L，IOA均值为0.60。相较于3DHED模拟结果，8个代表性站点IOA平均值提升了39.5%，RMSE平均值降低了6.9%，表明数据同化方法能够有效提高模型的模拟精度。

图3 蓝藻生物量观测值、3DHED模拟结果和DA模拟结果对比

3.3 改进数据同化结果分析

图4为8个湖区代表性站点蓝藻生物量观测值、DA模拟结果以及mDA模拟结果的对比。从图4可以看出，DA模拟结果和mDA模拟结果均能较为准确地模拟蓝藻生物量的变化趋势，但mDA捕捉到蓝藻生物量峰值的次数有一定增加，表明采用DA方法进行数据同化时，部分同化时刻遥感数据存在较大误差，mDA方法可进行高质量遥感数据的筛选，有效提升模拟精度。mDA模拟结果的RMSE值介于0.79～1.68 mg/L，8个站点RMSE均值为 1.21 mg/L，在DA模拟结果基础上降低了9.7%，mDA模拟结果的IOA值介于0.53～0.89，均值为0.68，在DA模拟结果基础上提高了13.3%。

图4 蓝藻生物量观测值、DA模拟结果和mDA模拟结果对比

图5为太湖30个采样点3DHED模拟结果、DA模拟结果和mDA模拟结果的IOA值及RMSE值对比。3种模拟方式下30个站点的IOA均值分别为0.43、0.64和0.71，RMSE均值为1.42 mg/L、1.27 mg/L和1.16 mg/L。各站点IOA值由大到小顺序为mDA、DA、3DHED，RMSE值由小到大顺序为mDA、DA、3DHED。总体而言，mDA的蓝藻生物量模拟精度最高，但少部分站点存在mDA模拟结果的RMSE高于DA模拟结果的现象。

(a) IOA值

表3为各个湖区3DHED模拟结果、DA模拟结果以及mDA模拟结果的IOA值和RMSE值的对比情况。从表3可以看出，DA模拟结果的IOA均高于3DHED模拟结果的IOA，而mDA模拟结果的IOA值在DA模拟结果的基础上又得到了进一步提升。对于3种模拟结果的RMSE，mDA模拟结果的RMSE值最低，模拟效果最好，整体上太湖各湖区RMSE值表现为西北部湖区大于东南部湖区。

4 讨论

在真实的环境下，与蓝藻相关的一些物理参数(如蓝藻沉降率和被捕食率等)是难以精确模拟和确定的，并且参数会随着时间而变化，其过程远比模型模拟复杂得多[29]。本文采用3DHED模型模拟2014—2016年太湖全湖蓝藻生物量的时空变化时，虽然能够在一定程度上反映太湖蓝藻的变化趋势，但总体模拟效果欠佳，峰值模拟效果有待进一步提高(图3)。蓝藻从生长到形成水华会经历下沉—休眠—复苏—上浮4个阶段[30]，冬季沉于底泥中的藻类细胞在春季适宜条件下迅速生长，于夏季聚集上浮形成水华[31]。水动力条件也是决定水华暴发的关键因子[32]，当蓝藻生物量聚集到一定程度后，风浪的扰动使蓝藻细胞通过碰撞形成大群体，风浪作用一旦减弱，便上浮形成表面可见水华[33]。结合本文模拟结果，3DHED模拟过程中未能考虑冬季沉底部分的蓝藻，同时对于风力作用下蓝藻的迁移过程描述不够清晰，这可能是导致模型对于蓝藻生物量峰值捕捉不够理想的重要原因[7]。采用遥感数据对3DHED进行EnKF数据同化后，由于各个观测时刻的蓝藻生物量得到及时更新，使得DA模式下的蓝藻生物量模拟精度得到明显提升(图3)，其IOA均值较3DHED模拟结果增加了48.8%，RMSE均值降低了10.4%。

表3 8个湖区3DHED、DA和mDA模拟结果的IOA值和RMSE值的对比

模拟误差和观测误差是数据同化的关键因素[16]，对同化性能有重要影响。设置较小的观测误差时，观测数据精度更高，同化结果会更接近观测值，反之，同化结果会更接近模拟值。一般情况下，观测数据的误差更小，更接近于真实值。李渊等[34-35]在进行不同同化试验时采用了较小的观测误差，取得了较为理想的试验结果。但遥感数据受到云覆盖等因素的影响，在某些时刻、某些区域的数据误差较大，在进行数据同化时，由于模型采用的仍是较小的观测误差，会导致遥感数据误差被低估，同化结果趋近于遥感数据，反而会更加偏离真实值。本文采用改进的数据同化模式后，结果显示mDA的蓝藻生物量模拟效果较DA整体上有了一定的提升，对蓝藻生物量峰值的捕捉频次也有所增加(图4)，全湖IOA均值在DA基础上增加10.9%，RMSE均值降低了8.6%，说明未改进前遥感数据在某些时刻存在误差被低估的现象，模型模拟结果趋向于误差被低估的遥感数据，从而偏离真实值。

太湖蓝藻水华总体呈现东南低西北高的趋势，这与入湖河流污染负荷输入密切相关。西南部的天目山水系，城镇较少，人口密度低，水质较好。西北部的南溪水系水流多来自工农业发达的平原和城镇地区，入湖污染负荷较大[36]。同时太湖夏季盛行东南风，受到太湖风浪的影响，蓝藻水华易在西北部聚集和堆积[37]，因此，重度水华主要集中在太湖西北部沿岸(如竺山湖、梅梁湾等)[38]。本文采用mDA方法进行太湖蓝藻生物量的模拟时，除了在一定程度上减小遥感数据误差对同化结果的影响外，还能较好地模拟蓝藻生物量的空间分布。结合图1和表3可以看出，RMSE值总体表现为西北部湖区大于东南部湖区，这是由于蓝藻生物量越大的区域，其RMSE值也越大。

湖泊藻类动态模型采用常微分方程描述湖泊生态系统的食物链动态和营养物质循环导致的自身结构不确定性[25]，使得模拟误差会随着模拟时间的增加而不断积累。数据同化技术通过实时加入观测数据，可及时更新初始场数据[39]，将一个初始场的使用时间缩短到一个观测周期内，有效降低了模拟误差的积累，此时初始场数据的误差对模型模拟精度起到决定性作用。遥感数据具有时间上较为连续、监测范围广等特点[40]，但若在观测时刻受到云雾等因素影响，提供的初始场数据误差可能大于模型模拟得到的初始场数据。本文以两种模式下提供的初始场数据模拟到下一观测时刻，比较得到的模拟结果与实测数据的RMSE，再筛选出与真实情况误差较小的初始场数据进行模型模拟，从而进一步提高模拟精度。近年来，国内学者针对太湖蓝藻生物量模拟开展了一系列研究。张成成[41]利用三维水生态模型模拟了蓝藻生物量的变化，验证站点的IOA值介于0.23～0.50。Huang等[28]采用EnKF方法通过更新模型状态变量和模型参数以提高太湖富营养化模型对浮游植物的模拟精度，其最佳模拟效果下IOA值达到了0.67。本研究mDA模式下的蓝藻生物量模拟结果改善效果明显(图5)，全湖平均IOA值为0.71，表明本文提出的mDA同化方法具有较高的模拟精度，尤其对于监测站点稀疏的大型湖库水华的精确模拟预测具有广阔的应用前景。值得注意的是，存在少部分站点mDA结果精度低于DA的现象(图5中拖山、洑东等站点)，这主要是由于本研究是通过计算整个湖区36个站点处同化结果与实测数据的RMSE值，来判断是否加入遥感数据进行同化，少部分站点遥感数据质量较差，但在整个湖区范围内精度较高，使得部分站点模拟效果不佳。因此后续的研究可以通过划分不同湖区来进行遥感数据分区的筛选，从而进一步提高数据同化精度。同时，太湖面积广阔，不同湖区的蓝藻水华严重程度亦有差别，进行不同湖区模型参数的分区率定对模型精度提高具有一定的效果。