基于GPU和数据同化的深水湖库水温与溶解氧中短期预报

孙博闻，宗庆志，杨晰淯，张袁宁，高学平

(天津大学 水利工程仿真与安全国家重点实验室，天津 300350)

1 研究背景

深水湖库通常存在季节性温度分层现象，水温影响着水体中生化反应[1-5]，由温度分层引起溶解氧分层与季节性缺氧还会诱发水环境水生态问题[6-8]，进而威胁供水安全及下游生态系统健康。在气候变化与水库调控的共同影响下，上述变化更加趋于复杂[9-10]，因此研究水温和溶解氧的变化对湖库生态环境安全十分必要。

当前对水库水温和溶解氧的研究大多关注其长期演变规律[11-13]，但在未来1～10天的中短期时间尺度上预报湖库水温与溶解氧变化对管理者同样具有重要意义[14]。如Weber等[15]利用水温-深度变化数据动态确定Grosse Dhuenn水库取水深度，在优化下泄水温同时还抑制了水库缺氧的发展；强风驱动产生的剧烈涌升流使日内瓦湖部分湖区的上层水温骤降，直接影响湖区溶解氧分布与下游供水[16]；余晓等[17]对潘家口水库全年的水温和溶解氧浓度进行时空分布特征和演化成因的研究，得出水温分层控制溶解氧分层的结论；刘畅等[18]发现大流量调度过程引起的强对流条件也可以产生短期削弱甚至完全消除由水库分层导致的缺氧现象。可见对湖库水温与溶解氧在中短期时间尺度进行预报，能够在调控湖库下泄水温这类常规管理中发挥作用，还能在应对极端天气等威胁供水与生态安全的应急管理工作中提供决策支撑。

数学模型是研究深水水库水动力水质演化规律的重要方法，但受制于模型计算速度以及模型初始条件与参数取值等原因，模型在中短期时间尺度上的预报效率与精度还有较大的提升空间。在提升模型计算效率方面，并行计算是当前主流解决方案。如能够采用合适的优化方法利用CPU(Central Processing Unit)处理数学模型逻辑运算部分，利用GPU(Graphic Processing Unit)执行并行运算部分，有可能大幅缩减模型运行时间。但利用GPU运算通常需要进行大量的代码改写，这也是造成当前大多数数学模型并未广泛采用GPU并行计算的重要原因。在提升模型计算精度方面，有学者提出利用数据同化的方法将观测数据与模型计算结果相融合，结合现场观测和数值模拟的优点，实现模型状态和参数不断更新，以有效提高模型的模拟精度[19]。但受限于初始条件、边界条件等不确定性的影响，目前数据同化方法在深水水库中短期水质预报的应用依然较少。

本文以引滦入津工程—大黑汀水库为研究对象，利用水库长时段、高频次自动监测数据为同化数据源，采用集合卡尔曼滤波算法作为同化方法，基于CE-QUAL-W2模型(the two-dimensional，laterally averaged，hydrodynamic and water quality model)建立水库水动力水质数学模型，选用GPU并行计算方法提高模型计算效率，构建大黑汀水库水温与溶解氧的数据同化系统，以期实现对水库垂向水温与溶解氧的动态变化的高精度、高效率的中短期预报，为水库水环境管理提供可靠的技术支撑。

2 深水湖库温-氧中短期预报计算

2.1 研究区域大黑汀水库位于唐山市迁西县城北5 km的滦河干流，地处北温带，四季分明，控制流域面积35 100 km2，总库容3.37亿m3。大黑汀水库是引滦入津枢纽工程的源头水库，承担着向唐山及天津供水的重要任务。自1990年代以来，大黑汀水库开始出现大规模的网箱养殖活动(图1)，人工投入的养殖饲料、鱼类排泄物等外源有机污染持续被投入，对水质造成了严重污染。水库管理局自2016年开始大面积清理养殖网箱，污染负荷被清除。但几十年养殖活动产生的污染物依旧沉积在水库底泥当中。对水库水温和溶解氧的监测数据表明，库区底部存在明显的缺氧区，水环境质量明显下降。

图1 大黑汀水库概况图Fig.1 General situation of Daheiting Reservoir

2.2 CE-QUAL-W2模型本文采用CE-QUAL-W2构建大黑汀水库水动力水质数学模型。CE-QUAL-W2是一个考虑纵向和垂向的二维水动力水质模型(简称W2模型)，模型假设水体横向平均，因此适宜模拟河道型水库的水动力与水质变化。模型边界条件包括上游流量、水温及水质、下游水位等，气象条件数据包括气温、露点温度、风速、风向和云量等。模型建立和参数率定过程参考姚嘉伟研究成果[20]，本文不再赘述。

2.3 基于OpenACC的GPU并行计算OpenACC(Open Accelerators)是一种新型的GPU并行方式，它利用导语向编译器提供编译和优化代码的相关信息，由于其执行效率与实现方法综合较优，近年来开始应用于水动力数学模型的GPU并行计算[21-23]。在对程序进行并行计算时，由于外部逻辑循环包括前后处理及串行代码很难实现在GPU上的并行计算，所以选择适当的循环代码进行并行计算有利于提升程序并行加速效果。其次，选出循环复杂的代码部分，在代码的开始和结尾处添加OpenACC导语代码，对程序进行编译，查看编译器的反馈信息并记录初步的计算用时。此时有可能会发生计算用时变长的情况，这是因为开始GPU并行计算后，增加了GPU和CPU之间的数据传输的耗时，如果传输的数据量过大，就会增加计算耗时，进而影响程序的计算性能。本研究基于OpenACC采用CPU-GPU协同计算结构实现W2模型的并行计算(图2(a))，通过OpenACC计算导语将W2模型水动力计算中较复杂的部分分配给GPU处理器，CPU处理器完成剩余计算任务及逻辑控制，选用Tesla K40c处理器作为GPU加速设备。由于数据同化过程需要执行大量模型计算，因此通过比较不同计算导语和模型网格数量对GPU加速效果的影响，确定数据同化阶段采用的GPU并行计算方案。

2.4 基于OpenDA的数据同化与预报OpenDA(Open data assimilation)是一个可快速实现参数校准和滤波测试等功能的数据同化工具箱。它提供了一个允许算法和模型进行简单交换的平台，通过将数据同化算法作为单独组件，降低数据同化的应用成本，近年来广泛应用于水文模型与水动力模型的数据同化领域[24-27]。本研究基于OpenDA进行水温与溶解氧的数据同化和预报，流程如图2(b)所示。其中数据同化系统包括预测模型系统与观测系统，同化算法选用集合卡尔曼滤波算法(EnKF，Ensemble Kalman filter algorithm)。在同化阶段预测模型都要运行集合数次，当有观测数据时，通过EnKF算法进行数据同化。数据同化后的结果被称为分析值集合，该集合的均值被认为是最优估计值，分析值集合将作为下个时间步预报的初始值，直到下一个有观测数据的时刻进行数据同化。

图2 计算流程Fig.2 Computing process

3 结果与分析

3.1 GPU加速效果

3.1.1 不同导语选择对GPU优化效果影响 表1列出了两种导语组合格式，相应计算用时对比如图3(a)所示。可以看出在GPU并行计算后，由于增加了主机和GPU之间的数据传输的耗时，因此不同导语组合方式的计算时间有较明显的差异。在不同计算网格数下，导语组合二的计算耗时始终小于导语组合一。随着网格数的增加，两组导语组合的计算耗时差逐渐增大，网格数为1000时耗时相差最大，达到了83 min。因为导语组合二在导语组合一的基础上，增加了copy子语和copyin子语，其中copy子语确保了数据由主机到GPU之间的双向传输的可行性，而copyin子语则可以将目标计算数据复制到显卡内存中，并且在GPU完成第一部分并行计算后不把数据传输回主机，所以可以有效的减少数据传输消耗时间，提升GPU并行效率。因此选择适当的OpenACC导语，可以有效减少数据传输所消耗的时间，提升GPU并行加速效果。

表1 OpenACC的不同导语组合

3.1.2 不同网格数量下GPU加速效果分析 本研究网格数设定为300、350、400、450、500、550、600、650、700、750、800、850、900、950和1000。选取3.1.1节中GPU并行效率较高的导语组合二作为本实验的OpenACC导语形式，实验结果如图3(b)所示。在计算网格数小于600时，CPU单精度的计算模式耗时最少，但随着网格数的增加，GPU并行计算耗时低短的优势开始逐渐显现出来。而当计算网格数增加到950时，虽然仍存在数据传输的耗时问题，但由于GPU强大的并行计算能力，CPU-GPU协同计算比CPU串行计算节省用时可达46 min。为保证GPU加速效果，同时考虑模拟精度要求，最终确定模型纵向划分为69个单元(Segment)，垂向最多33层(Layer)。

图3 不同导语组合或网格数量条件下的计算用时对比Fig.3 Comparison of computational time under different combinations of directives or grid quantity

3.2 水温及溶解氧数据同化

3.2.1 同化和预报性能判别标准 采用大黑汀水库坝前断面水质监测站点2019年实测5个月的温度和溶解氧数据(采样间隔30 min)，对CE-QUAL-W2模型的状态变量和参数同时进行同化。以均方根误差(RMSE)、一致性系数(IOA)和百分比偏差系数(PBIAS)作为评价标准(式(1)—式(3))，RMSE越接近0表示模拟效果越佳，IOA用来表征模拟值和观测值的一致性，取值在0～1之间，越接近1表示二者之间的一致程度越高，PBIAS结果以百分比定量化给出模拟值比观测值整体被低估或高估的平均趋势，PBIAS为0表示模拟效果最佳，PBIAS>0表示模拟值倾向于低估，反之PBIAS<0表示模型倾向于高估。

(1)

(2)

(3)

3.2.2 同化效果分析 图4为有无数据同化的情况下，坝前断面水下1 m、7 m和13 m处的模拟水温与实测水温、模拟溶解氧和实测溶解氧间的偏差值。可以看出，在经过EnKF方法对模型参数和状态变量同时进行更新后，模型的模拟精度有了很大的提升，数据同化模型在坝前断面水下1 m、7 m和13 m处的水温平均模拟偏差分别为1.01 ℃、0.31 ℃和0.21 ℃，较无数据同化的模拟结果分别提升了68.8%、51.6%和41.2%；溶解氧平均模拟偏差分别为0.38 mg/L、0.22 mg/L和0.18 mg/L，较无数据同化的模拟结果分别提升了71.2%、63.6%和54.4%。

图4 W2模型模拟结果与数据同化系统模拟结果相比实测水温和溶解氧的偏差值Fig.4 The water temperature errors and the dissolved oxygen (DO) errors

3.3 水温及溶解氧中短期预报

3.3.1 水温中短期预报 以数据同化系统为基础，考虑边界条件等不确定性前提下，运行W2模型进行水库水温预报。参考Zhang等[28]的研究结果，以单日内是否所有剖面数据都监测出温跃层为依据判断分层状态，即5月27日—10月4日为分层稳定期，此前为分层形成期，此后为分层削弱期。为分析数据同化系统的水温预报性能，绘制分层期、混合期与全时段下各水深处未来1～10天的水温预报结果，如图5(第1—3列)所示。可以看出，在1 m与13 m处分层期的预报性能优于混合期，7 m处两阶段预报性能各有优劣。从RMSE指标来看，随着预报时段由1天延长至10天，预报误差由0.22～0.35 ℃增大至0.77～1.09 ℃，即预报精度会随着预报时段延长逐渐降低。其他学者的研究结果也表明，预报时长与预报精度显著相关[29-31]，预报时长的增加会使预报结果的不确定性增大。从IOA指标来看，分层期7天预报期内各深度的IOA值均超过0.65，具有较好的一致性；预报期延长至10天后，7 m处的IOA由0.78降至0.57，降幅远超1 m与13 m深度。从PBIAS指标来看，1 m和13 m处预报均倾向于高估，7 m处则倾向于低估，但高/低估均在±3%以内。近年来引滦工程常采用潘家口水库大流量调度抑制大黑汀水库缺氧发展[19-20，32]，水库调度提高了水温预报的不确定性，在温跃层附近的水温预报精度依然存在提升空间。从各项统计指标结果可以看出，构建的数据同化系统能够较好地预报水库1～10天内水温变化。

图5 分层期、混合期与全时段各水深处未来1～10天的水温预报结果评价Fig.5 Evaluation of water temperature forecast results for the next 1-10 days at different depths

3.3.2 溶解氧中短期预报 为分析数据同化系统对溶解氧的预报性能，在水库分层状态和混合状态下对坝前断面各水深处未来1～10天的溶解氧预报结果进行评价，结果如图6。随着预报时段由1天延长至10天，RMSE由0.06～0.25 mg/L增大至0.45～1.09 mg/L，可见随着预报时段延长，不确定因素的累计效应对预报会产生一定影响。分层期7天预报期内各深度的IOA值均超过0.65，具有较好的一致性；预报期延长至10天后，13 m处的IOA由0.8降至0.64，降幅远超1 m与7 m深度。从PBIAS指标来看，1 m和13 m处预报均倾向于高估，7 m处则倾向于低估，但高/低估均在±3%以内。从各项统计指标结果可以看出，构建的数据同化系统能够较好的预报水库1～10天内溶解氧变化。

图6 各水深处未来1～10天的溶解氧预报结果评价Fig.6 Evaluation of DO forecast results for the next 1-10 days at different depths

4 结论

本文采用集合卡尔曼滤波算法作为同化方法，基于CE-QUAL-W2模型建立水库水动力水质数学模型，基于OpenACC的GPU并行方法提升模型计算效率，构建大黑汀水库水温与溶解氧的数据同化系统，在中短期时间尺度上实现对水库水温与溶解氧动态变化的高精度、高效率预报，主要结论包括：(1)GPU并行计算中，不同导语组合的计算耗时差异明显，网格数量为1000时，计算耗时相差可达83 min。当OpenACC采用kernels导语加copy和copyin子语的组合时，GPU并行加速的表现最好。当网格数超过600后，网格数量越多，GPU并行计算的效率优势越大。(2)采用EnKF方法进行数据同化，同时同化状态变量与模型参数，不同深度处水温模拟结果较无同化条件提升了41.2%～68.8%，溶解氧模拟结果提升了54.4%～71.2%。有助于提升对水温和溶解氧的中短期预报性能。(3)1 m、7 m和13 m处水温预报的平均偏差分别在0.60 ℃、0.32 ℃和0.28 ℃以内；1 m、7 m和13 m处溶解氧预报平均偏差分别在0.54 mg/L、0.32 mg/L和0.24 mg/L以内。预报结果符合水库水温与溶解氧的中短期变化规律，能够为大黑汀水库的供水与生态安全提供支撑。