一种基于通用计算的流域水质实时模拟方法

陈 军，李 婷，禹政阳

(成都信息工程大学资源与环境学院，成都 610225)

0 引 言

水质模型模拟污染物在水环境中的运动和迁移，是水生态环境评价、水污染预警预报的重要工具。近几十年来，国内外许多学者在水质模拟方面已开展了大量的研究工作[1-11]。二维水质模型研究污染物在水平方向上的扩散，常用于大型河流的水质模拟，逐渐成为流域水质模拟的重点研究方向。

二维水质模型将流域按一定的空间规则划分为大量不相互重叠的空间网格，通过网格的污染物迁移模拟污染物的空间扩散；将污染物扩散过程按一定时间间隔分解为大量的时间片元，通过时间片的迭代模拟污染物的时间变化规律。因此，二维水质模型属于计算密集性运算，随着网格划分空间精度和时间片元的较小，运算量呈倍增加。传统的CPU计算难以实现高精度的二维水质模拟。由于GPU强大的并行计算能力，将通用计算引入到水质模拟成为未来的趋势[12-15]。然而，在目前显卡硬件条件下，一次送入GPU的数据大小有一定限制。对于大型流域，网格空间尺度精细到一定尺度后，GPU很难一次性读入。因此，如何实现流域范围的水质实时模拟，成为水质模拟必须解决的关键问题。

另一方面，流域的前期水情是影响水质模拟精度的重要因素。由于数据获取能力的限制，流域的前期水动力和水质数据往往难以获取。为解决该问题，目前大多数模型采用“预热”的方式，即利用有限的观察数据对模型进行多年的迭代，将迭代结果作为流域前期水情，典型的如SWAT[16-20]，这种方式最大问题是性能问题；另一部分模型要求输入河道所有位置上的水深和污染物浓度数据，典型的如EFDC[21-23]。由于难以获取河道所有位置上的水深和浓度，一般采用空间离散的观测点插值得到。空间插值的误差则成为该方法中前期水情误差的主要来源。

为此，本文提出一种基于通用计算的流域水质实时模拟方法。通过降低模拟网格规模和水质模拟的并行化设计，以提升水质模拟性能；通过引入稳定流场[24-28]，以提升水质模拟的精度。实验证明，相对于传统水质模拟，其性能提升显著，模拟精度达到要求，其成果具有一定的理论意义和实践价值。

1 通用计算框架下的流域网格划分

1.1 栅格模式的流域网格划分

在二维水质模型中，基本计算单元称为“微元”，它代表了水体的一个微小区域。“微元”之间无缝连接，以模拟研究区域的水体。流域的水质模拟采用微元的污染物模拟计算实现。为便于通用计算，本文将“微元”设定为规则的矩形网格，即对流域按相等的分辨率划分为一定数量的行和列构成的矩阵，矩阵中每一个单元为二维水质模型中的“微元”，也是栅格模式下的一个像元。

栅格模式的划分虽然便于通用计算处理，但随着分辨率的提高，像元数量的增加，存储和计算量成倍增加。如果将流域的整个矩形区域的像元均放置于显存，在现有技术下所能表达的流域范围和分辨率有限。

1.2 模拟网格的定义

为河道设置最大水体范围，并假定在模拟时段内水和污染物运动局限于该范围内的像元。以此为基础，将流域范围内的像元分为两类，一类为水体像元，另一类为外部像元。由于现实中水体像元仅占全流域范围的一小部分，如果仅将水体像元放入显存中运算，则在相同的显存下支持更大范围和更高分辨率的流域水质模拟。

然而，在二维水质模型中，每一个像元均需要与周围8邻域像元进行水和污染物交换，需要将与水体邻接的外部像元的相关属性同时存储到显存中。将与水体邻接的外部像元称为边界像元，则参与水质模拟的网格包括水体像元和边界像元两部分，如图1所示。其中，水体像元为实际水质计算的网格单元，边界像元为邻接的水体像元提供邻接属性的网格单元。

图1 流域水质模拟网格划分Fig.1 The classification of cells for watershed water quality simulation

模拟网格集合M表示为：

M={celln| 1≤n≤Nm}

(1)

式中：Nm为集合中模拟网格的个数；n为各模拟网格在集合中存储的序号，用于唯一标识每一个模拟网格。

在通用计算中，水质模拟所需的地形、水深、水速、污染物浓度均表达成长度为Nm的数组，模拟网格的数量决定了流域的水质模拟所需的显存大小。表1给出了实验流域不同空间分辨率下模拟网格数量。从表1可见，模拟网格数量不及网格总数的1%，并且随分辨率的增加，占比不断减小，并无限接近水体占整体流域的比例。因此，在流域中定义模拟网格，降低了水质模拟计算的显存需求，为流域尺度的水质模拟通用计算奠定了基础。

表1 不同栅格分辨率条件下的模拟网格数量及占比Tab.1 The number and proportion of simulated grids under different grid resolution conditions

1.3 网格邻域拓扑的构建

模拟网格在污染物扩散与迁移计算时，需要周围8邻域网格的属性信息。因此，在模拟网格划分的基础上，还要构建网格邻域拓扑。对于每一个模拟网格，其8邻域按照图2所示的顺序编码，依次记录邻域网格在模拟网格集合中的序号。

图2 中心网格的8邻域编码Fig.2 8 Neighborhood coding of the center cell

为便于水质模拟的通用计算，模拟网格的邻域拓扑也采用数组表示。由于一个模拟网格有8个邻域编码，邻域拓扑数组的大小是Nm的8倍。为便于邻域网格的检索，邻域拓扑数组按模拟网格序号的顺序依次存储。

2 流域水质模型的并行化设计

2.1 通用计算对水质模型的基本要求

可并行性是通用计算的基本要求。在通用计算中，一个任务被划分为大量的子任务，这些子任务的执行过程具有相对独立性。为实现基于通用计算的水质模拟，水质模拟过程也应该进行可并行性分解。水质模拟包括时间迭代和空间迭代。时间迭代是从上一个时间片的污染物状态得到下一个时间片的模拟结果，时间片之间是串行的。因此，水质模型的并行化关键在于空间迭代。栅格像元之间的计算独立性，是水质模型通用计算的关键。而栅格像元之间的计算独立性，要求在模型中，模拟网格仅能修改它自身的属性，不能修改其他网格的属性。

2.2 模拟网格的属性化

为实现水质模拟的通用计算，将模拟网格的属性cell表示为一个多元组：

cell={H,s,b,Wf,vf,cf,Wt,vt,ct}

(2)

式中：H为网格的地形高度；s为网格类型，分为水体像元和边界像元两类；b为网格的邻域拓扑信息；Wf、vf、cf分别表示时间片上网格的初始水深、水速和污染物浓度；Wt、vt、ct分别表示时间片上网格的终止水深、水速和污染物浓度。

在时间片上，依据模拟网格的{Wf,vf,cf}，计算得到该时间片的{Wt,vt,ct,}。下一个时间片计算时，将其进行交换。这种表示方法确保了时间片之间的数据独立性。

2.3 模拟网格的水动力并行算法设计

水动力模型是水质模型的基础。通过水动力模型，模拟出网格的水深和水速，依据水深和水速，模拟污染物的扩散、迁移和衰减。在一个时间片上，模拟网格的水动力模拟核心在于水深、水速的计算。对于流域任一模拟网格，在模拟时表示为中心网格r，其8邻域网格表示为b。时间片上的水动力模拟按如下方式进行。

2.3.1 模拟网格水深的计算

在栅格模式下，模拟网格的水体看作为一个立方体水柱。由于模拟网格为栅格模式下的一个像元。模拟网格的水量V由式(3)计算：

V=WC2

(3)

式中：W为模拟网格的水深；C为像元分辨率。

(4)

网格供水主要途径为降水和支流入汇。在仅考虑干流河道网格的水质模型中，上游干流入流、支流入汇、排污口均被假设为一个位于河道的水源点。由于流域中降水最终通过上游来水、支流入汇方式进入模拟河段，而河道上的降水量对于整个流域而言仅占很小一部分，因此模型中忽略降水项。对于河道的水源点，以流量为单位输入模型。设某水源点流量为ξ，m3/s，则时间片上水源点所在网格因供水增加的水深ξr为：

(5)

式中：N为一个小时的时间片数。

(6)

式中：(Δx、Δy)为邻域网格坐标相对中心网格坐标的偏移，像素；vmax表示时间片上允许的最大速度标量。

依据式(6)，当模拟网格水速等于vmax时，网格水流100%流出。因此，水质模拟的时间片间隔Δt可采用下式来计算[29]：

(7)

2.3.2 模拟网格的水速计算

(8)

式中：α为与地表摩擦力、重力加速度和水体密度等相关的常量。

(9)

2.3.3 时间片上水动力模拟的并行计算流程

为保证水动力计算的可并行性，时间片上的水动力模拟分两步进行。第一步，水速增量计算。首先依据水源点供水量计算模拟网格的水量增量，然后进行水速增量模拟，计算各网格单元新的水深和水速；第二步，水流迁移计算。以第一步得到的网格单元新水深和水速为基础，进行水流迁移计算，得到的水深和水速作为下一个时间片的初始水深和水速。

2.4 模拟网格的水质模拟并行算法设计

将水质模拟与水动力模型耦合，在水速增量、水流迁移环节中加入污染物迁移、扩散和衰减计算。

2.4.1 水速增量环节中的水质模拟

(10)

在水速增量环节中，同时考虑污染物降解。污染物降解系数以天为单位，需要将其换算到一个时间片上污染物的降解率。设污染物降解系数在流域内为一个常数kd，则时间上污染物降解率kr为：

kr=1-(1-kd)1/(24×N)

(11)

2.4.2 水流迁移环节中的水质模拟

(12)

2.4.3 污染物扩散模拟

在二维水质模型中，中心网格r向邻域网格b扩散，导致中心网格污染物减少；邻域网格b向中心网格r扩散，导致中心网格污染物增加。将污染物扩散系数设置为一个常量E，则相邻两网格单元因扩散导致的污染物实际迁移量与它们的浓度差有关。

将相邻两网格单元的水均看作为一个立方体水柱，水柱的高度等于水深，并认为污染物扩散发生在两水柱交叠区域，其水柱高度Wbr为：

Wbr=max[min(Hr+Wr,Hb+Wb)-max(Hr,Hb),0]

(13)

式中：Hr、Hb分别表示中心网格和邻域网格的地形高度，m。

图3 污染物扩散交叠区Fig.3 The contaminant diffusion overlap zone

为定量描述水流混合模拟导致的污染物浓度变化，将E换算到时间片上污染物向外扩散的距离e：

(14)

(15)

式中：W为流域出口控制点水深；下标表示模拟小时数；n为当前已模拟的小时数；h为前推小时数；A为累积误差阈值。

为提高模型模拟精度，利用实际水位和浓度数据对模型参数进行率定，率定后vmax=2 m/s、α=0.12、σ=0.995、kd=0.05、E=0.08。将牡丹江中游2012年5月1日0时各水源点瞬时流量代入模型中，基于流域90 m地形数据创建稳定流场。各水源点瞬时流量和浓度如表2所示。

表2 不同栅格分辨率条件下的流域水质模拟时间Tab.2 The time of water quality simulation at different grid resolutions

从表2可见，模拟网格总数在20 000以内时，GPU和CPU均可完成流域水质模拟，GPU比CPU快10倍左右；随着流域分辨率提高，网格规模及每小时迭代次数成倍增加，CPU耗时越来越多，仅能通过几小时水质模拟所需时间估算2年水质模拟总时间，而GPU仍能满足流域水质实时模拟要求。

3 基于稳定流场的流域水质模拟

3.1 稳定流场的创建

为合理模拟流域前期水情，引入稳定流场以解决由空间离散的观测点水深和浓度插值生成初始边界条件造成的误差问题。假定流域存在若干个水源点，其流量恒定，通过径流汇流模型不断迭代，直至流域出口处的水深、水速和浓度最终趋于稳定。显然，现实中不存在稳定水源，稳定流场也不可能存在。但在一个具体的时间点上，流域干流河道上的干流入流点、支流入汇点和排污点均存在瞬时流量和浓度。如果以模拟时段的起始时刻各水源点的瞬时流量作为稳定供水量，瞬时浓度作为稳定污染源排放浓度，通过模型迭代创建的稳定流场可近似作为模拟时段的前期水情。在稳定流场创建过程中，采用式(16)定量描述稳定流场是否创建完成。

(16)

式中：W为流域出口控制点水深；下标表示模拟小时数；n为当前已模拟的小时数；h为前推小时数；A为累积误差阈值。

为提高模型模拟精度，利用实际水位和浓度数据对模型参数进行率定，率定后vmax=2 m/s、α=0.12、σ=0.995、kd=0.05、E=0.08。将牡丹江中游2012年5月1日0时各水源点瞬时流量代入模型中，基于流域90 m地形数据创建稳定流场。各水源点瞬时流量和浓度如表3所示。

表3 起始时间各水源点瞬时流量和化学需氧量浓度Tab.3 The instantaneous flow rate and COD concentration of each water source at the start time

在h=10 h、A=0.05的条件下，叠加表4所示的水源点初始流量和浓度，迭代200 h后稳定流场生成，如图4所示。

图4 实验流域稳定流场创建过程Fig.4 The stable flow field creation process of experimental basin

图5为牡丹江中游2012年5月1日0时稳定流场下模拟的初始化学需氧量浓度分布图。将水质监控断面的模拟初始浓度、IDW插值的初始浓度与实测浓度对比，结果如表4所示。从表4可见，稳定流场模拟的流域初始浓度分布更贴近实际情况，为流域高精度水质模拟奠定了基础。

表4 水质监控断面初始浓度模拟精度对比Tab.4 The comparison of simulation accuracy of initial concentration of water quality monitoring sections

3.2 水质模拟精度评价

将牡丹江中游2012年5月1日至2014年4月30日干流入流点、支流入汇点和排污点逐小时的流量和浓度代入稳定流场进行模型迭代，得到模拟时段内逐小时的河道化学需氧量浓度分布。采集该流域温春大桥和海浪水质监测断面栅格单元逐小时的模拟浓度，制作浓度变化曲线，与实测浓度和EFDC模型(已进行最优化参数率定)分别进行对比，如图6所示。

图6 水质监测断面模拟浓度与实测浓度对比Fig.6 The comparison of simulated concentration and measured concentration of water quality monitoring sections

从图6可见，本文模型对温春大桥断面和海浪河断面的浓度变化趋势拟合效果均达到EFDC模型相近的效果。将断面的实测浓度cr与模拟浓度cm按下式计算拟合度R：

(17)

式中：N为断面实测浓度序列个数。

断面模拟精度如表5所示。从表5中可见，两个监测断面的模拟精度均达到85%以上，与EFDC模型具备相近的模拟精度。

表5 模拟时段内水质监测断面模拟精度Tab.5 The simulation accuracy of water quality monitoring sections during simulation period

4 结 语

由于二维水质模拟属于计算密集型运算，高精度网格划分和计算性能的瓶颈之间难以调和的矛盾，迫切需要引入通用计算来提升流域水质模拟的性能和效率。本文从水质模拟的性能和模拟精度两方面进行了探讨。在水质模型性能方面，以通用计算为框架，针对流域空间尺度，将河道及边缘栅格作为模拟网格，降低了网格规模，为流域精细化水质模拟奠定了基础。之后，针对通用计算对数据并行化要求，将水质模拟的各个阶段，包括水速增量、水流迁移和污染物扩散模拟，均进行了并行化设计。实验发现，将通用计算引入到水质模拟，可极大地提高模拟性能，使流域水质实时模拟成为现实。

为提高水质模拟精度，引入了稳定流场。在水质模拟时，首先将初始时刻流域水源点的流量和浓度进行迭代，直至流域出口的水深和浓度稳定为止，将模拟出来的稳定流场作为水质模拟的前期水情进行模拟时段的水质模拟。对比发现，本文的前期水情模拟方法相对于传统的空间插值具有更高的可信度。

最后，通过牡丹江中游2年的水质模拟，证明了在模拟网格相当的条件下，本文方法将水质模拟时间从传统水质模拟的4.5 h缩短到几分钟，并且还达到了传统水质模拟相近的精度，这为流域水质实时模拟奠定了基础。

□