CN109583035A - 基于元胞自动机的城市面源污染物累积冲刷过程计算方法 - Google Patents

Abstract

本发明公开了一种基于元胞自动机的城市面源污染物累积冲刷过程计算方法，包含：构建元胞地理场景，获得每个元胞的下垫面类型、高程和随时间变化的降雨量；根据下垫面下渗参数，进行产流计算和流向计算，进而通过汇流计算获得每个元胞的当前水深；根据污染物累积和冲刷相关系数和干期天数，计算每个元胞的污染物累计量，并进行污染物输移过程的模拟计算，获得每个元胞的陆相污染物转移量和水相污染物输移量，进而计算水相污染物浓度。本发明结合水利学和污染物累计冲刷原理，具有明确的物理含义和理论基础，能够有效模拟污染物的累计冲刷过程，为城市面源污染累积冲刷模拟提供了新方法，为城市面源污染负荷削减和控制提供指导和依据。

Description

基于元胞自动机的城市面源污染物累积冲刷过程计算方法

技术领域

本发明涉及环境科学领域，尤其涉及一种基于元胞自动机的城市面源污染物累积冲刷过程计算方法。

背景技术

随着城市化的快速发展，城市不透水面不断扩大，降雨径流冲刷汇流所形成的面源污染已成为城市水环境质量下降的主要原因。面源污染模型能够对污染物的迁移转化过程进行定量模拟，已成为城市面源负荷来源解析和治理的重要研究内容。目前城市面源污染模型在模拟污染物冲刷时大多采用黑箱方法或灰箱方法，只能获得研究区的整体冲刷量而不能反映局部冲刷特征，模型参数缺乏物理含义，对污染物冲刷过程的机理描述尚不充分。

元胞自动机是一种时间、空间、状态都离散的网格动力学模型，具有模拟复杂系统时空演化过程的能力，尤其适合变量存在局部非连续或者剧烈梯度变化的情景，在城市复杂地形面源污染累积与冲刷过程模拟中具有独特的优势。

因此，充分考虑污染物冲刷时的降雨溅蚀作用、径流冲刷作用和颗粒物再沉积等物理作用机制的基础上，结合水力学和污染物累积冲刷原理，构建城市面源污染物累积冲刷的元胞自动机模型，实现基于元胞自动机的城市降雨径流和污染物冲刷过程时空动态模拟，对于城市面源污染负荷削减和控制具有重要的指导意义，这也本发明的主要研究课题。

发明内容

发明目的：本发明目的是提供一种基于元胞自动机的城市面源污染物累积冲刷过程计算方法，以更机理的描述污染物的累积冲刷过程，为面源污染累积冲刷过程模拟提供新方法。

技术方案：本发明所述的基于元胞自动机的城市面源污染物累积冲刷过程计算方法，主要包括如下步骤：

(1)构建元胞地理场景模型，包括定义地理场景模型元胞大小、领域关系、模型边界和出口，根据获取的研究区的下垫面类型、DEM高程和随时间变化的降雨量，获得每个元胞的下垫面类型、高程和随时间变化的降雨量；

(2)根据设置的下垫面下渗参数利用初损后损法进行每个元胞的产流和下渗量的计算得到元胞的净雨深，并根据中心元胞与邻域元胞高程的水位差进行流向计算，进而通过汇流计算获得每个元胞的当前水深；

(3)根据设置的污染物累积和冲刷相关系数和干期天数，采用指数累积公式计算每个元胞的干期污染物累积量，并考虑污染物的降雨溅蚀、径流冲刷和再沉积作用，对污染物输移过程进行模拟计算，得到径流冲刷强度、降雨冲刷强度和污染物再沉降强度，从而获得每个元胞的陆相污染物转移量和水相污染物输移量，进而计算水相污染物浓度。

进一步的，考虑到模型计算的结果有效性，所述方法还包括(4)结合实测的流量和污染物数据，利用纳什效率系数(EF)和皮尔逊相关系数(CC)进行计算结果的评价。

进一步的，所述步骤(2)中利用初损后损法将下渗简化为初损和后损两个部分，降雨开始到出现产流时的水量损失为初损量(I₀)；产流后的水量损失为后损，一般用后损平均入渗率(f)来表示；通过设置初损量和平均后损入渗率计算净雨深(h)，也即每个元胞的水深，计算公式如下：

h＝X-I₀-ft_c-X_n

式中，h为X形成的元胞的水深(净雨深)，mm；X为总降雨量，mm；I₀为初损量，mm；f为平均后损入渗率，mm/min；t_c为超渗产流历时，min；X_n为后损阶段非超渗历时的降雨量，mm。

进一步的，所述步骤(2)中方法进行水流流向的判断，具体为：比较中心元胞与邻域元胞高程的水位差，排除高于平均高程的邻域元胞，循环该过程直到没有高于平均高程的邻域元胞可以被排除，确定可以进行水量转移的邻域元胞。

进一步的，所述步骤(2)中汇流计算具体为：在产流和流向确定的基础上，采用曼宁公式计算一个时间步长内中心元胞可以转移的水量，同时根据水位差得出可以转移的潜在最大水量，将二者相比，得出一个时间步长内中心元胞转移的真实水量。计算公式如下：

式中，T为中心元胞水深转移至邻居元胞所需时间，既水量/污染物转移时间步长；V为流速，m/s；h为中心元胞水深，m；s为中心元胞与第i个邻居元胞间的水面坡度；n为曼宁糙率系数，s/m^1/3；t为模型时间步长；D为中心元胞与邻居元胞i的距离，m，对于四邻域而言，D为元胞边长，对于斜向邻域元胞，D为倍的元胞边长；f_i为水深，表示中心元胞向邻居元胞i转移的水量。

比较时间步长t与计算出的T的大小，计算中心元胞向邻居元胞i转移的水量，如果T小于时间步长t，认为水量f_i可以全部转移；否则只能转移相应比例的水量。根据元胞平均高程与邻居元胞的高差计算最多可以转移至的水量，每个邻居元胞转移水量f_i求和可以得出转移总水深。并将转移总水深与中心元胞水深进行对比，确定每个元胞的实际转移水量。将上述规则应用于整个元胞空间，计算出每一时刻整个元胞空间内所有元胞的转移水量和剩余水量，并根据下一时刻步长的降雨和下渗量，获得下一时刻末的元胞水量，循环计算，直到整个模拟过程结束。

进一步的，所述步骤(3)中干期污染物累积量计算公式如下：

式中，M为干期沉积层污染物累积量，单位面积或者边沿长度的质量；td为干期时长，天；C₁为最大增长可能，单位面积或者边沿长度的质量，C₂为增长速率常数，1/天。

进一步的，所述步骤(3)中按照如下公式进行污染物冲刷过程计算，得到每个元胞陆相污染物转移量：

d_i＝v_sic_i

其中h₀＝0.33D_R

Ω_e＝F(Ω-Ω₀)，其中Ω＝ρ_wgS_fq

Δp_i＝(r_ri+e_ri-d_i)*T*L²

式中，m_T为单位面积的沉积层颗粒物总质量，kg/m²；m_i为单位面积的沉积层中第i组颗粒物质量，kg/m²；P为降雨强度，m/s；ρ_s和ρ_w分别为沉积物密度和水密度，kg/m³；g为重力加速度常数，m/s²；v_si为第i组颗粒物在水中沉降速度，m/s；c_i为第i组颗粒物在水体的浓度，kg/m³；a_d为沉积层颗粒物可分离性，kg/m³；a_d0为沉积物初始可分离性，kg/m³；q为单宽流量，m²/s；h为元胞水深，m；h₀为临界水深，m；b为常数；D_R为雨滴平均直径，mm；Ω_e为有效径流功率，W/m²；Ω为计算的总径流功率，W/m²；Ω₀为临界径流功率，W/m²；S_f为摩阻比降；F为比例系数；r_ri为径流冲刷强度；e_ri为降雨冲刷强度；d_i为污染物再沉降强度；△p_i为每个元胞陆相污染物转移量；T为水量转移时间步长，也即模型污染物模拟时间步长，s；L为元胞边长，m。

进一步的，所述步骤(3)中根据每个元胞初始污染物累计量M和陆相污染物转移量△p_i，利用初始污染物累积量减去陆相污染物转移量即可得出此时刻每个元胞陆相污染物总量，初始水相污染物量加上陆相污染物转移量即可得到水相污染物总量。根据每个元胞的水相污染物量除以元胞总的水量，即可以求出每个元胞的水体污染物浓度。在水量转移的基础上，按照每个元胞此时刻的污染物浓度，进行污染物量转移计算，得出整个元胞空间每个元胞当前时刻末的污染物量和污染物浓度，循环计算，直到整个模拟过程结束。

有益效果：与现有技术相比，本发明综合考虑了降雨溅蚀作用、径流冲刷作用和颗粒物再沉积作用等物理机制，建立了一种面源污染物累积冲刷过程元胞自动机模型。该模型考虑了颗粒物粒径对冲刷的影响，使用了真实的时间步长，相较于传统模型，本研究提出的模型具有明确物理含义和理论基础，不仅能够准确模拟面源污染物的累积冲刷输移过程，而且能够模拟不同级配颗粒物的冲刷规律和集水区出口的颗粒物粒径组成。

附图说明

图1是本发明的一个实施例的流程示意图；

图2是本发明的元胞边界条件示意图；

图3是流向判断算法流程示意图；

图4是污染物冲刷计算流程示意图；

图5是流量验证结果示意图；

图6是总悬浮物TSS验证结果示意图。

具体实施方式

如图1所示，本发明实施例公开的一种基于元胞自动机的城市面源污染物累积冲刷过程，包括以下步骤：

(1)获取研究区的下垫面类型、数字高程数据(DEM)和随时间变化的降雨量，设置元胞大小、时间步长、模型边界和出口，建立元胞基础地理场景。

该步骤具体包括：

(1-1)获取研究区的下垫面类型和数字高程数据(DEM)构建元胞基础地理场景，获得每个元胞的下垫面类型(例如：道路和屋顶等)和高程，设置模型的边界和模型的出口，定义每个元胞的大小、曼宁糙率系数、坐标、邻域关系，模型的时间步长和模拟时长，相关参数如表1所示。如图2所示，本发明中元胞的邻域关系为Moore型八邻域，边界元胞设置为极大值从而保证水量无法进入，本发明暂不考虑城市管网的影响和蒸散发损失。

表1地理元胞模型参数集示例表

参数	意义	说明
			CA_length	元胞边长	模型设置
n	曼宁糙率	下垫面材质平均糙率
			T_stop	模拟时长	根据降雨过程进行设置
t	时间步长	设置模拟时间的步长

(1-2)利用雨量计测定或通过气象站点获得研究区的降雨量数据，利用Arcgis软件进行站点降雨量数据的时空插值，获得每个元胞不同时刻的降雨量数据。

(2)根据设置的下垫面下渗参数(初损量、后损平均入渗率等)，进行产流计算和流向计算，进而通过汇流计算获得每个元胞的当前水深。

该步骤具体包括：

(2-1)产流计算：利用初损后损法将下渗简化为初损和后损两个部分，降雨开始到出现产流时的水量损失为初损量(I₀)；产流后的水量损失为后损，一般用后损平均入渗率(f)来表示；通过设置初损量和平均后损入渗率计算净雨深(h)，也即每个元胞的水深，计算公式如下：

h＝X-I₀-ft_c-X_n

式中，h为X形成的元胞的水深，也即净雨深，mm；X为总降雨量，mm；I₀为初损量，mm；f为平均后损入渗率，mm/min；t_c为超渗产流历时，min；X_n为后损阶段非超渗历时的降雨量，mm。

比较降雨量和初损量，当累积降雨量达到初损量时，元胞开始具备产流条件。根据模型时间步长，分别计算出当前时间步长内各个元胞的降雨量和后损入渗量，并比较大小。若降雨量大于入渗量，该元胞产流，且产流量为降雨量与入渗量的差值；反之，该元胞不产流。将此规则用于整个元胞空间，结合元胞的大小，即可得到所有元胞当前时刻的水深。

(2-2)流向判断：根据DEM和元胞水深进行水流流向的判断，确定中心元胞与邻域元胞的平均高程，比较邻域元胞与元胞平均高程，排除高于平均高程的邻域元胞，循环该过程，直到没有高于平均高程的元胞可以被排除，确定可以进行水量转移的元胞，流向判断算法流程图如图3所示。

(2-3)汇流计算：在产流和流向确定的基础上，采用曼宁公式计算一个时间步长内中心元胞可以转移的水量，同时根据水位差，得出可以转移的潜在最大水量，将二者相比，可得出一个时间步长内中心元胞转移的真实水量。

具体的，假设中心元胞的水共流入k个邻居元胞，用第i个表示在k个邻居元胞中的某一个。首先，根据曼宁公式计算中心元胞水深h转移至第i个邻居元胞所需时间T，(也既水量/污染物转移时间步长)，计算公式如下：

其中，V为流速，m/s；h为中心元胞水深，m；s为中心元胞与第i个邻居元胞间的水面坡度；n为曼宁糙率系数，s/m^1/3；D为中心元胞与邻居元胞i的距离，m，对于四邻域而言，D为元胞边长，对于斜向邻域元胞，D为倍的元胞边长。

进而，比较时间步长t与计算出的T的大小，计算中心元胞向邻居元胞i转移的水量(用水深f_i表示)。如果T小于时间步长t，认为水量可以全部转移，否则只能转移相应比例的水量。水深f_i使用如下公式计算：

根据元胞平均高程与邻居元胞i间的高差计算最多可以转移至元胞i的水量，如果高差大于f_i，则水量转移为f_i；若高差小于f_i，则转移水量的量使得两元胞高程相同，没有水位差。每个邻居元胞转移水量f_i求和可以得出转移总水深(f₀)，将其与中心元胞水深进行对比，若中心元胞水深h小于f₀，即中心元胞没有足够的水进行转移，每个邻居元胞水量转移量需按照水深进行缩小，重新计算f₀；若中心元胞水深h大于f₀，转移水量不变。将此规则应用于整个元胞空间，计算出每一时刻整个元胞空间内所有元胞的转移水量和剩余水量。

通过以上计算可以得到当前时刻末每个元胞的水量，当前时刻末每个元胞的水量作为下一时刻初始时元胞的水量，根据下一时刻步长的降雨和下渗量，获得下一时刻末的元胞水量，循环计算，直到整个模拟过程结束。

(3)根据设置的污染物累积和冲刷相关系数和干期天数，计算每个元胞的污染物累计量，并进行污染物输移过程的模拟计算，获得每个元胞的陆相污染物转移量、水相污染物输移量和水相污染物浓度。

该步骤具体包括：

(3-1)污染物累积过程计算：设置污染物累积干期天数，干期天数可以通过气象局网站公布的每日降水量数据进行确定。采用应用较多的SWMM软件中的指数累积公式计算每个元胞的污染物累积量；其中，干期污染物累积量计算公式如下：

式中，M为沉积层污染物累积量，单位面积或者边沿长度的质量；t为干期时长，天；C₁为最大增长可能，单位面积或者边沿长度的质量，C₂为增长速率常数，1/天。C₁和C₂为经验系数，可通过文献获得。

(3-2)污染物冲刷过程计算：结合元胞自动机的转换规则进行污染物输移过程的模拟计算，计算径流冲刷强度r_ri、降雨冲刷强度e_ri和污染物再沉降强度d_i，进而计算每个元胞陆相污染物转移量(△p_i)，具体计算公式如下：

d_i＝v_sic_i

Δp_i＝(r_ri+e_ri-d_i)*T*L²

其中，

Ω_e＝F(Ω-Ω₀)，Ω＝ρ_wgS_fq

式中，m_T为单位面积的沉积层颗粒物总质量，kg/m²，可由(3-1)中计算的污染物累积量进行计算；m_i为单位面积的沉积层中第i组颗粒物质量，kg/m²，可利用粒径分析仪测量粒径组分比例测定，并乘以总的污染物量M来获得；P为降雨强度，m/s；ρ_s和ρ_w分别为沉积物密度和水密度，kg/m³，为经验参数；g为重力加速度常数，m/s²；v_si为第i组颗粒物在水中沉降速度，m/s，可采用Stokes公式进行计算；c_i为第i组颗粒物在水体的浓度，kg/m³；a_d为沉积层颗粒物可分离性，kg/m³；a_d0为沉积物初始可分离性，kg/m³，为模型经验参数；q为单宽流量，m²/s；h为元胞当前时刻水深，m；h₀为临界水深，m，为模型经验参数；b为常数，单位为1，为模型经验参数；D_R为雨滴平均直径，mm；Ω_e为有效径流功率，W/m²；Ω为计算的总径流功率，W/m²；Ω₀为临界径流功率，W/m²，为模型经验参数；S_f为摩阻比降；F为一个比例系数，反映了有效径流功率与超过临界径流功率部分的总径流功率的比值；T为水量转移时间步长，也即模型污染物模拟时间步长，s；L为元胞边长，m。

根据计算获得每个元胞的陆相污染转移量△p_i，判断△p_i的正负，△p_i若为正，说明累积污染物从陆相被冲刷入水中量大于沉降量，污染物净进入水体；△p_i若为负，说明污染物沉降量大于被冲刷进入水中的量，污染物从水体输出。

根据每个元胞陆相污染物转移量△p_i，得出此时刻每个元胞水相污染物总量P_水i(初始水相污染物量加上陆相污染物转移量)和陆相污染物总量P_沉i(初始陆相污染物量减去陆相污染物转移量)；由水相污染物量除以元胞总的水量，可以求出每个元胞水体中污染物浓度。在水量转移的基础上，按照每个元胞此时刻的污染物浓度，进行污染物量转移计算，得出整个元胞空间每个元胞当前时刻末的污染物量和污染物浓度，循环计算，直到整个模拟过程结束，污染物冲刷计算过程如图4所示。

(4)结合实测的流量和污染物数据，利用纳什效率系数(EF)和皮尔逊相关系数(CC)进行计算结果的评价。

使用流量计进行降雨过程流量的测量，获得每个时刻模型出口的实测流量数据。污染物数据获取时，利用经过校正的聚乙烯塑料瓶在出水口收集全过程的径流，考虑到降雨的初期冲刷效应以及雨量的大小，前期水样收集的时间间隔应设置相对较短(例如：每1-2min收集一次)，后期水样收集时间间隔可适当加大，每次降雨过程一般收集15-20瓶水样，记录每瓶水样收集满所用时间。按照采样顺序进行编号并在进行样品的实验室分析。水质检测分析方法按照参照国家标准(《水和废水监测分析方法(第四版)》)进行，总悬浮物(TSS)检测采用重量法。

根据获得实测的流量和TSS数据，计算纳什效率系数(EF)和皮尔逊相关系数(CC)评价模拟结果的可靠性和适用性。其中，EF和CC为已有成熟的计算公式，不作具体介绍。流量验证结果如图5所示，污染物验证结果如图6所示。

以上所揭露的仅为本发明一种较佳实施例而已，不能以此来限定本发明之权利范围，因此依本发明权利要求所作的等同变化，仍属本发明所涵盖的范围。

Claims (8)

1.基于元胞自动机的城市面源污染物累积冲刷过程计算方法，其特征在于，该方法包括：

(1)构建元胞地理场景模型，包括定义地理场景模型元胞大小、领域关系、模型边界和出口，根据研究区的下垫面类型、DEM高程和随时间变化的降雨量，获得每个元胞的下垫面类型、高程和随时间变化的降雨量；

(2)根据设置的下垫面下渗参数利用初损后损法进行每个元胞的产流和下渗量的计算得到元胞的净雨深，并根据中心元胞与邻域元胞高程的水位差进行流向计算，进而通过汇流计算获得每个元胞的当前水深；

(3)根据设置的污染物累积和冲刷相关系数和干期天数，采用指数累积公式计算干期污染物累积量，并考虑污染物的降雨溅蚀、径流冲刷和再沉积作用，对污染物输移过程进行模拟计算，得到径流冲刷强度、降雨冲刷强度和污染物再沉降强度，从而获得每个元胞的陆相污染物转移量和水相污染物输移量，进而计算水相污染物浓度。

2.根据权利要求1所述的基于元胞自动机的城市面源污染物累积冲刷过程计算方法，其特征在于：所述步骤(2)中利用初损后损法将下渗简化为初损和后损两个部分，降雨开始到出现产流时的水量损失为初损量I₀；产流后的水量损失为后损，用后损平均入渗率f表示；每个元胞的净雨深计算公式如下：

h＝X-I₀-ft_c-X_n

式中，h为X形成的元胞的水深，即净雨深；X为总降雨量；t_c为超渗产流历时；X_n为后损阶段非超渗历时的降雨量。

3.根据权利要求1所述的基于元胞自动机的城市面源污染物累积冲刷过程计算方法，其特征在于：所述步骤(2)中水流流向判断方法，具体为：比较中心元胞与邻域元胞高程的水位差，排除高于平均高程的邻域元胞，循环该过程直到没有高于平均高程的邻域元胞可以被排除，确定可以进行水量转移的邻域元胞。

4.根据权利要求1所述的基于元胞自动机的城市面源污染物累积冲刷过程计算方法，其特征在于：所述步骤(2)中汇流计算具体为：在产流和流向确定的基础上，采用曼宁公式计算一个时间步长内中心元胞可以转移的水量，同时根据水位差得出可以转移的潜在最大水量，将二者相比，得出一个时间步长内中心元胞转移的真实水量；计算公式如下：

式中，T为中心元胞水深转移至邻居元胞所需时间，既水量/污染物转移时间步长；h为中心元胞水深；s为中心元胞与第i个邻居元胞间的水面坡度；n为曼宁糙率系数；t为模型时间步长；D为中心元胞与邻居元胞i的距离，对于四邻域而言，D为元胞边长，对于斜向邻域元胞，D为倍的元胞边长；f_i表示中心元胞向邻居元胞i转移的水量；

每个邻居元胞转移水量f_i求和可以得出转移总水深，将转移总水深与中心元胞水深进行对比，确定每个元胞的实际转移水量；计算出每一时刻整个元胞空间内所有元胞的转移水量和剩余水量，并根据下一时刻步长的降雨和下渗量，获得下一时刻末的元胞水量，循环计算，直到整个模拟过程结束。

5.根据权利要求1所述的基于元胞自动机的城市面源污染物累积冲刷过程计算方法，其特征在于：所述步骤(3)中干期污染物累积量计算公式如下：

式中，M为干期沉积层污染物累积量，单位面积或者边沿长度的质量；td为干期时长；C₁为最大增长可能，单位面积或者边沿长度的质量，C₂为增长速率常数。

6.根据权利要求1所述的基于元胞自动机的城市面源污染物累积冲刷过程计算方法，其特征在于：所述步骤(3)中按照如下公式进行污染物冲刷过程计算，得到每个元胞陆相污染物转移量：

d_i＝v_sic_i

其中h₀＝0.33D_R

Ω_e＝F(Ω-Ω₀)，其中Ω＝ρ_wgS_fq

Δp_i＝(r_ri+e_ri-d_i)*T*L²

式中，m_T为单位面积的沉积层颗粒物总质量；m_i为单位面积的沉积层中第i组颗粒物质量；P为降雨强度；ρ_s和ρ_w分别为沉积物密度和水密度；g为重力加速度常数；v_si为第i组颗粒物在水中沉降速度；c_i为第i组颗粒物在水体的浓度；a_d为沉积层颗粒物可分离性；a_d0为沉积物初始可分离性；q为单宽流量；h为元胞水深；h₀为临界水深；b为常数；D_R为雨滴平均直径；Ω_e为有效径流功率；Ω为计算的总径流功率；Ω₀为临界径流功率；S_f为摩阻比降；F为比例系数；r_ri为径流冲刷强度；e_ri为降雨冲刷强度；d_i为污染物再沉降强度；△p_i为每个元胞陆相污染物转移量；T为水量转移时间步长，也即模型污染物模拟时间步长；L为元胞边长。

7.根据权利要求1所述的基于元胞自动机的城市面源污染物累积冲刷过程计算方法，其特征在于：所述步骤(3)中根据每个元胞初始污染物累计量和陆相污染物转移量，利用初始污染物累积量减去陆相污染物转移量得出每个元胞陆相污染物总量，初始水相污染物量加上陆相污染物转移量得到水相污染物总量；根据每个元胞的水相污染物量除以元胞总的水量，求出每个元胞的水体污染物浓度；在水量转移的基础上，按照每个元胞的污染物浓度，进行污染物量转移计算，得出整个元胞空间每个元胞的污染物量和污染物浓度，循环计算，直到整个模拟过程结束。

8.根据权利要求1所述的基于元胞自动机的城市面源污染物累积冲刷过程计算方法，其特征在于：所述方法还包括：步骤(4)结合实测的流量和污染物数据，利用纳什效率系数和皮尔逊相关系数进行模拟计算结果的评价。