基于链式反应的湖泊型流域非点源污染多级靶向治理目标确
定方法
技术领域
[0001] 本
发明涉及
水利工程领域，尤其涉及水环境管理领域，具体为基于链式反应的湖泊型流 域非点源污染多级靶向治理目标确定方法。
背景技术
[0002] 近年来，随着点源污染逐渐得到控制，非点源污染成为影响水环境的主要污染源。大量 的泥沙、氮磷污染负荷进入江河、湖库，会引起
水体悬浮物浓度升高、有毒有害物质增加、 溶解
氧减少、水体出现富营养化及
酸化趋势，这不仅会打破
生态系统原有的平衡状态，还会 影响人类正常的生产生活，威胁人体健康。从目前来看，非点源污染产生的氮磷营养物质已 成为世界范围内普遍存在的重大水环境问题，并成为当今水环境保护的主要障碍因子。识别 和量化面源污染特征是合理保护流域水环境的科学
基础，这不仅符合当前的学术研究热点， 而且在实际应用中也具有重要的理论指导意义。
[0003] 有研究表明，非点源污染控制工作能否取得长足效果，取决于能够准确
定位和系统修复 流域非点源污染靶向治理目标。流域非点源污染靶向治理目标(targeted treatment objectives， TTOs)又称优先控制区(priority management areas,PMAs)或者(critical source areas,CSAs)， 被定义为虽然拥有相对较小的流域面积占比，但是对目标水域贡献着不成比例的高污染负荷 区域，它是高污染源区、高产流区以及高
土壤侵蚀区的综合表征。从本质上来看，非点源污 染的靶向治理目标等同于点源污染中的重大污染源，点源污染中的重大污染源指排污量大且 对环境影响剧烈的重点工矿企业，但相比于点源污染的排放
位置相对固定、排放规模可准确 评估且排放
频率受人为管控等特征，非点源污染具有产生的不确定性、排放的随机性、污染 源难以追踪且影响因素众多。以全流域尺度开展非点源污染治理工作既不科学，也不经济。 因此，在开展非点源污染治理工作中，不是绝对的进行全流域治理才能确定取得显著的水环 境改善效果，相反，若能将有限的人
力和财力资源放置在面积虽小但是对目标水质能产生极 大影响的靶向治理目标区域内，则会取得事半功倍的效果。总体上来看，靶向治理目标的确 定是为了寻求最大的水质改善效果与最小的经济投入之间的平衡点。
[0004] 在我国，非点源污染的定量化计算是流域污染物总量控制的重要内容。但是就目前来看， 所实施的污染物目标总量控制
框架是
自上而下、指标分解、工程落实的方式实现负荷的削减， 缺乏建立流域污染物
排放量与目标水体水质之间的响应关系，使得污染物总量削减与水体水 环境
质量目标之间处于脱节的状态，这一
瓶颈势必会对我国当前和未来水环境精细化管理， 由总量管理向质量管理转型发展造成一定负面影响。
[0005] 相比于河道型流域而言，湖泊型流域拥有更大的水域面积，对入湖污染负荷拥有更大的 自净能力。然而，从目前的研究成果来看，多是以子流域的污染负荷产生量，或者以陆域与 水体结合部之间的污染负荷输出量作为对湖泊水环境影响的评判标准。对于湖泊型流域而言， 受纳水体包括入湖河网水系和湖泊水体两部分，只有全面系统的认识河网水系对于非点源污 染输移转化的影响和湖泊水体水质对于污染物输入的响应，才能实现湖泊型流域非点源污染 的精细化管理。
发明内容
[0006] 本发明的目的是提供一种基于链式反应的湖泊型流域非点源污染多级靶向治理目标确定 方法。本发明的目的是通过以下方案实现的：
[0007] 本方法为一种
基于链式反应的湖泊型流域非点源污染多级靶向治理目标确定方法，该方 法针对现有的流域非点源污染靶向治理目标识别过程中存在的问题，提出将目标水体水质的 响应作为流域非点源污染靶向治理目标识别的基准，将污染物的产生、河网水系输移影响过 程以及水体水质响应过程看成一个完整的水环境系统来看待，能够克服传统的流域非点源污 染模拟过程中单一获取总出口污染负荷量以及流域非点源污染靶向治理目标识别过程中缺乏 对河网水系和目标水体水质考虑的不足，能够为高效的制定流域非点源污染治理方法提供坚 实的理论指导作用。该方法需要收集流域气象、下垫面情况以及水文水质监测数据，具体包 括以下步骤：
[0008] 1)构建流域分布式水文、非点源污染模型。以流域地理高程、土地利用类型、土壤类型 以及坡度变化作为基本输入数据，建立以水文响应单元(Hydrology Respond Units，HRUs) 为基本计算网格，以子流域为统计单元的流域分布式水文、非点源污染模型，采用重要控制 断面、子流域出口以及流域总出口的水文、水质监测数据进行模型率定工作，获取适合于目 标流域水文、营养物循环的关键参数，该部分内容是后续计算工作的基础。
[0009] 2)构建子流域空间拓扑关系矩阵R。河网水系是流域水文循环的重要组成结构，也是非 点源污染物质迁移转化的重要途径。因此，在流域水文和非点源污染模拟过程中，各子流域 之间都是以河网水系产生水力联系的。即：①相邻的子流域之间具有明确的上下游隶属关系； ②水量和污染物质以子流域作为传输媒介，具备清晰的理论基础。对于相邻的两个子流域Wi和Wj，Wi为紧邻Wj的上游流域，其水量和污染物质由Wi转移到Wj的条件概率为1，反之， 条件概率为0，
[0010] 3)构建子流域污染输出负荷矩阵E。根据已建立的流域分布式水文、非点源污染模型， 模拟计算子流域内的污染物质的产生、迁移和转化情况。子流域是由多个HRU单元组成， HRU单元是由相似的土地利用类型、土壤和坡度组成，某一子流域Wi的污染物输出量Ei是 该子流域范围内所有HRU单元污染物输出量的总和。
[0011] 4)构建河网水系输移影响关系矩阵F。相比于流域尺度来看，河道的面积相比较小，但 是河道对于污染物的迁移转化起到了很重要的作用，如河道
沉积物对于污染物的
吸附，藻类 对于污染物的分解、转化。考虑到污染物在河道内的物理、化学和
生物过程较为复杂。根据 各子流域的污染负荷输入及输出结果，采用污染物通量相对变化值来表征某一段河网水系对 于污染物输移的影响，河网水系输移影响系数的取值范围在0到1之间，影响系数的数值越 小，说明河网水系对于污染物输移影响的程度越大，反之，河网水系对于污染物输移影响的 程度越小。
[0012] 5)计算湖泊水体水质污染物响应系数k。相比于河道型流域而言，湖泊具有更大范围的 水体面积和更强的自净能力，如污染物质的物理沉降、污染物之间的化学反应以及浮游生物 的吸收与释放等。污染物质进入湖泊水体后，很快混合均匀，可以将湖泊水体视为一个混合 均匀的反应器，并且不考虑渗漏，湖泊内污染物质仅发生衰减反应并且符合一级反应动力学。 对于一个长期存在的湖泊，其水量、污染物循环已处于一个平衡稳定状态。为了表征湖泊水 体对于污染物质的一系列复杂的物理、化学和生物过程，采用水体一维完全混合降解模型获 取湖泊水体水质对于污染物质的响应系数k。
[0013] 6)构建子流域入湖污染贡献量矩阵D。基于上述建立的子流域空间拓扑关系矩阵R、子 流域污染输出负荷矩阵E、河网水系输移影响关系矩阵F以及湖泊水体水质污染物响应系数 k，计算各子流域入湖污染贡献量Di。对于具有n个子流域的入湖河流而言，污染物最多经过 n次转移后达到入湖口，再通过大水体的稀释、降解作用后，形成影响湖泊水环境状态变化 的污染贡献量。
[0014] 7)确定湖泊水环境评估周期T。降水是流域水文循环的主要驱动力，地表产汇流和湖泊 水量是污染物质迁移转化的载体。受大气环流的影响，降水在时间序列上表现出一定的周期 性。对于一个稳定的湖泊流域而言，降水的周期性变化势必会
对流域的产水、产污带来很大 程度上的影响，因此，湖泊水环境状态的变化也存在着周期性。采用水文小波分析方法确定 湖泊水环境评估周期T，这是进行流域非点源污染多级靶向治理目标识别的重要前提和基础。
[0015] 8)构建目标湖泊水质达标率P——湖区水体污染物超标通量α——流域污染负荷累积量 的映射β关系。基于灰色概率分布，建立湖泊水环境状态变化、湖区水体污染物超标通量变 化和流域污染负荷变化一一映射关系。为了定量化评价湖区水质的变化情况，
结合水功能区 目标水质，引入目标水质达标保证率P的概念，即在某一评估周期内，针对某一目标水质浓 度，水质达标时间出现的概率。与之相对应的，为了达到某一水质目标，需要对湖泊过量的 污染物通量进行削减，某一目标水质达标保证率P会相对应着湖区水体污染物超标通量α(t)。 进一步的，结合湖泊水体水质污染物响应系数以及流域入湖污染总负荷量，计算湖区水体污 染物超标通量所对应的流域污染负荷累积比例β(％)。同时，为了更贴近流域非点源污染治 理工作，引入流域非点源治理措施实施效率系数M，用以表征流域非点源污染潜在削减水平。
[0016] 9)识别不同水质达标率下的流域非点源污染多级靶向治理目标。根据步骤5)得到的各 子流域污染贡献量Di，为了消除子流域面积对于计算结果的影响，采用污染贡献量强度Ii对 各子流域进行排序，并绘制子流域污染贡献量累积曲线CD和子流域面积累积曲线CA。采用 步骤8)得到的不同湖泊水质达标率下所对应的流域污染负荷累积比例作为流域非点源污染 靶向治理目标识别的划分依据，获取为了达到某一湖泊水质达标率时需要优先控制的子流域 位置、面积等信息。
[0017] 进一步的，步骤1)中，选取流域分布式水文模型SWAT(Soil&Water Assessment Tool) 作为模拟流域水文、污染物循环的基础工具。其水量平衡、泥沙输移、氮磷循环的关键数学 公式如下：
[0018] ①SWAT模型水量平衡方程：
[0019]
[0020] 式中：SWt为土壤最终
含水量(mm)；SW0为土壤前期含水量(mm)；t为时间步长(day)； Rday为第i天降雨量；Qsurf为第i天的
地表径流(mm)；Ea为第i天的
蒸发量(mm)；Wseep为第i天存在于土壤坡面底层的渗透量和侧流量(mm)；Qgw为第i天
地下水出流量(mm)。
[0021] ②泥沙输移方程：
[0022] sed＝11.8(Qsurf·qpeak·areahru)0.56·Kusle·Cusle·Pusle·LSusle·CFRG (2)[0023] 式中，sed为土壤侵蚀量，t；Qsurf为地表径流量，mm/hm2；qpeak为峰值流量，m3/s；areahru为水文响应单元(HRU)的面积，hm2；Kusle为土壤侵蚀因子；Cusle为植被
覆盖和管理因子； Pusle为保持措施因子；LSusle为地形因子；CFRG为粗糙碎屑因子。
[0025]
[0026] 式中，Nnit,ly为第ly层土层中 转化为 的氮含量(kg N/ha)；frnit,ly为硝化过程中产生的 氮损失分数；frvol,ly为挥发过程中产生的氮损失分数；Nnit|vol,ly为第ly层土层硝化过程和挥发过 程共同引起的氮损失量(kg N/ha)。
[0028]
[0029] 式中，Ndeniy,ly为反硝化过程所产生的氮损失量(kg N/ha)；NO3ly为第ly层土层中的硝酸盐含 量(kg N/ha)；βdenit为反硝化速率；γtmp,ly为第ly层土层中营养物循环
温度因子；γsw,ly为 第ly层土层中营养物
循环水分因子；orgCly为土壤层中有机
碳含量(％)；γsw,thr为营养物循 环中发生反硝化过程的水分因子
阈值。
[0030] ④残留物有机磷的矿化：
[0031] Pminf,ly＝0.8·δntr,ly·orgPfrsh,ly (5)
[0032] 式中，Pminf,ly为为残留物有机磷经过矿化作用转化的无机磷(kg N/ha)，δntr,ly为残留物降解系 数，orgPfrsh,ly为第ly层土层中的残留物有机磷含量(kg/ha)。
[0033] 溶解态无机磷与活性态有机磷之间的磷平衡方程：
[0034]
[0035] 式中，Psol|act,ly为溶解态无机磷与活性态无机磷之间转化的磷含量(kg P/ha)；Psolution,ly为溶解 态无机磷含量(kg P/ha)；Pact,ly为活性态无机磷含量(kg P/ha)；pai为磷可利用指数。当Psol|act,ly大于0时，磷从溶解态无机磷向活性态无机磷转化；当Psol|act,ly小于
0时，磷从活性态无机磷 向溶解态无机磷转化。
[0036] SWAT模型中关于流域水文、营养物循环的详细机理过程描述参见以下著作 (S.L.NEITSCH,J.G.ARNOLD,J.R.KINIRY,J.R,et al.WILLIAMS.SOIL AND WATER ASSESSMENT TOOL THEORETICAL DOCUMENTATION.2009.Grassland Soil and Water Research Laboratory——Agricultural Research Service Blackland Research Center——Texas AgriLife Research)。
[0037] 进一步的，步骤2)中，定义一：Wi为对应子流域i，i为子流域编号，n为子流域总个 数，令{W(n),n≥0}为随机离散的子流域序列，I＝{i0，i1，i2，…}为离散的状态集合，若对 任意n个子流域和任意的i0，…，in-1，i，j∈I，满足条件概率：
[0038] P(Xn+1＝j|X0＝i0,···,Xn-1＝in-1,Xn＝i)＝P(Xn+1＝j|Xn＝i)＝pij (7)[0039] 则称这n个随机子流域序列{W(n),n≥0}满足
马尔科夫链过程，其本质含义是子流域i的状态 仅依赖与之紧邻的子流域j的状态，而与其余子流域i0，…，in-1的状态无关。
[0040] 定义二：子流域i和子流域j是紧邻的子流域，水量和污染物从子流域i转移到子流域j 的条件概率P(Wn+1＝j∣Wn＝i)＝pij称为空间马尔科夫链{W(n),n≥0}在子流域n的一步转移概 率，其中i，j∈I，n≥0，其矩阵形式为：
[0041]
[0042] 该转移概率矩阵具有如下的性质：
[0043] pij≥0 (9)
[0044]
[0045] 在实际应用中，需要根据流域的实际情况定义子流域空间拓扑矩阵、子流域污染物输出 矩阵、河道影响系数矩阵，具体的：
[0046] 对于n个子流域W1，W2，…，Wi，Wj，…，Wn，其空间拓扑矩阵可用n×n的矩阵进行 表示，Wi为紧邻Wj的上游流域，其水量和污染物质由Wi转移到Wj的条件概率为1，反之， 条件概率为0，其子流域空间拓扑关系矩阵数学表达式为：
[0047]
[0048]
[0049] 根据定义二，矩阵Rn×n的行向量之和为0或1。当行向量为1时，说明子流域有唯一确 定的出口；当行向量为0时，说明该子流域位于流域总出口处。在经过最多n次转移后，某 一子流域Wi的水量和污染物质将从流域总出口流出，即：
[0050]
[0051] 式中m为某一子流域Wi从源头到流域总出口所经历的最多子流域个数。
[0052] 进一步的，步骤3)中，子流域Wi的污染物输出量Ei是该子流域范围内所有HRU单元 污染物输出量的总和，其计算公式如下：
[0053]
[0054] 式中Ei为子流域Wi的污染物输出量，在SWAT模型中，Ei不仅考虑了陆域范围内随径流 进入水体的污染负荷，还包括了该子流域内河道部分对污染物的贡献量；AHRU(i,j)为子流域 Wi内HRUj的面积；CHRU(i,j)为子流域Wi内HRUj的污染负荷输出强度。
[0055] 根据上述获取的不同子流域污染物输出量Ei，定义并建立各子流域的污染物输出量矩阵， 其数学表达式为：
[0056] E＝(E1,E2,···,Ei,···,En)T (15)
[0057] 进一步的，步骤4)中，河道影响系数矩阵的数学表达式如下：
[0058]
[0059] 式中Fi为子流域i的河道影响系数，其取值范围在0到1之间，河道影响系数表示河道对 污染物的削减或者贡献作用，本质上是污染物在河道中迁移转化的效率。
[0060]
[0061]
[0062]
[0063] 式中Lout为子流域Wi出口断面的污染负荷量；Lin为子流域Wi入口断面的污染物负荷量， 其值为上游子流域Wi-1污染负荷的输出量Lout,Wi-1和子流域Wi的污染负荷输入量LWi之和。
[0064] 进一步的，步骤5)中，对于一个具有m个污染源(入湖河流)的湖泊型流域而言，湖 泊水体一维完全混合降解模型的数学表达式为：
[0065]
[0066] 式中c为湖泊污染物浓度(mg/L)；t为研究时段(d)；V为湖泊水量(m3)；qi为第i个 污染源(入湖河流)的入湖流量(m3/d)；ci为第i个污染源(入湖河流)入湖污染物浓度(mg/L)； Q为出湖水量(m3/d)，kj为某一污染物的综合降解系数(1/d)。
[0067] 对于一个长期存在且相对稳定的湖泊来说，若已知入湖水量、水质，出湖水量，湖泊水 量、水质的基础上，即可以求得某一污染物的综合降解系数，其数学表达式为：
[0068]
[0069] 在实际应用中，由于污染物的综合降解系数受到多种因素影响，如水量、水温等，因此需 要选取不同的典型水文年来计算污染物的综合降解系数。
[0070] 进一步的，步骤6)中，子流域入湖污染贡献量计算数学表达式为：
[0071] Di＝(Rn×n·Fn×n)i·Ei·kj (22)
[0072] 式中Di为某一子流域Wi的入湖污染物贡献量，kg，i≤n；Rn×n为子流域空间拓扑关系矩 阵；Fn×n为河道影响系数矩阵；Ei为某一子流域污染输出量，kg；kj为某一污染物的综合降 解系数，1/d。公式(22)的物理含义是子流域Wi，Wi+1，…，Wn-1，Wn的输出污染负荷经过 n-1，n-i-1，…，1，0次转移后到达入湖口，再经过大水体的稀释和降解作用后，最终形成影 响湖泊水环境状态的污染物贡献量。
[0073] 进一步的，步骤7)中，需要合理的划定湖泊水环境状态的评估周期，选取流域水文、 营养物循环的主要驱动力——降水作为评估周期划分的指标。采用小波分析方法分析降水量 的周期性，小波分析是一种窗口大小固定但形状可变的时频局部化分析方法。小波分析的关 键在于引入满足一定条件的基本小波函数Ψ(t)，指具有震荡特性、能够迅速衰减到零的一类 函数，经过伸缩和平移可得到一簇函数：
[0074]
[0075] 式中Ψa,b(t)为分析小波或连续小波，a为尺度伸缩因子，反映小波的周期长度；b为时间平移 因子，t反映时间上的平移。
[0076] 一般选取Morlet小波函数：
[0077]
[0078] 其中，c为常数，取6.2；i表示虚数。Morlet小波的伸缩尺度a与周期有如下关系：
[0079]
[0080] 对于上述给定的小波函数Ψ(t)，水文时间序列f(t)∈L2(R)的连续
小波变换为：
[0081]
[0082] 其中： 是Ψ(t)的复共轭函数；Wf(a,b)是小波变化系数。在实际研究中，时间序列常常 是离散的，如f(k·Δt)(k＝1,2,...,N；Δt为取样时间间隔)，则上式的离散形式为：
[0083]
[0084] 将时间域上的关于a的所有小波变换系数的平方进行积分，即为小波方差：
[0085]
[0086] 小波方差随尺度a的变化过程称小波方差图，反映了
波动的
能量随尺度的分布。通过小 波方差图，可以确定一个水文序列中存在的主要时间尺度，即主周期。
[0087] 水文周期小波分析方法的具体计算过程参见以下著作(魏凤英，现代
气候统计诊断与预 测技术，1999，气象出版社)。
[0088] 进一步的，步骤8)中，湖泊水环境状态变化、湖区水体污染物超标通量变化和流域污 染负荷变化一一映射关系结构如图10所示。
[0089] 水质达标保证率P的数学表达式为：
[0090]
[0091] 式中T达标为湖区水质达到目标水质的时间，T评估周期为评估周期时长，可为年、月、日。
[0092] 与之相对应的，湖泊污染物超标累积通量α(t)的数学表达式为：
[0093]
[0094] 式中C观测和C目标分别为湖区水质观测值(mg/L)和水质目标值(mg/L)，为了保险起见，以 年内月份最差的水质浓度作为该年湖泊污染物超标量的计算依据。V代表水位为典型代表水 文年下特征水位所对应的湖泊水量(m3)。
[0095] 与之相对应的，流域污染负荷累积比例β(％)的数学表达式为：
[0096]
[0097] 式中L总负荷量为典型代表水文年的流域入湖总污染负荷量(t)，其余符号含义如上所示。
[0098] 本发明方法区别于传统的流域非点源污染靶向目标识别方法，主要体现在两个方面： (1)建立子流域空间拓扑关系矩阵，考虑河网水系对于每个子流域污染物输出量的影响，计 算得到每个子流域的入湖污染贡献量，克服了以往传统方法单一获取流域总出口的入湖污染 负荷量，而忽略了各子流域的污染贡献量变化过程的不足；(2)构建起子流域污染输出量—— 河网水系影响——湖泊水体水质响应的完整驱动链条，真正将污染物质从产生、输出、迁移、 降解等一系列复杂的物理化学生物过程有机的结合起来，克服了以往传统方法忽略河网水系 的水力学联系以及水体水质的响应过程，而单纯的以子流域产生量或输出量作为流域非点源 污染评价基准的不足；(3)构建起湖泊水质达标保证率——湖区污染物累积通量——流域污 染物累积量的一一对应的映射关系，并根据不同的湖泊水质达标保证率识别流域非点源污染 治理的靶向目标，克服了以往传统方法进行流域非点源污染靶向治理目标识别时存在单一性、 “一刀切”的不足。该方法能够考虑目标水体水质对于子流域污染负荷的响应过程，将流域 非点源治理工作真正做到有的放矢。
[0099] 本发明的有益效果：
[0100] 在进行流域非点源污染靶向治理目标识别过程中，考虑了河网水系对于污染物输移的影 响以及大水体对于污染物的降解作用，识别的基准考虑了水体水质的响应变化过程，克服了 以往传统方法中主观的将污染产生量或污染输出量作为识别依据的不足。本方法提出的链式 反应结构，将污染物质产生、输出、迁移再到水体水质响应这一系列复杂的物理、化学和生 物过程有机的结合起来，将流域水环境看成一个完整的系统来对待。同时，从流域水环境精 细化管理的
角度来看，本方法通过构建流域污染物质的链式反应结构，可以精确识别出每个 子流域对于湖泊水环境状态的贡献量，扭转了以往粗放式的流域水环境管理方式。此外，本 方法可以识别出不同湖泊水环境状态下所对应的流域非点源靶向治理目标，流域相关部
门可 以结合当地经济发展条件和环境保护需求，实施多梯度的流域污染治理方案。综上所述，本 方法能够弥补已有的流域非点源污染治理工作中的不足，结合受纳水体水质的响应过程，精 准识别流域非点源污染优先治理区域，为实施污染削减最佳管理措施提供精确的靶向治理目 标。
附图说明
[0101] 图1为本发明的基于链式反应的湖泊型流域非点源污染多级靶向治理目标识别方法流程 图；
[0102] 图2为本发明的
实施例提供的洱海流域位置图；
[0103] 图3为本发明的实施例提供的洱海流域典型入湖河流水量模拟结果与实测数据率定期、 验证期对比图；
[0104] 图4为本发明的实施例提供的洱海流域典型入湖河流总氮模拟结果与实测数据率定期、 验证期对比图；
[0105] 图5为本发明的实施例提供的洱海流域典型入湖河流总磷模拟结果与实测数据率定期、 验证期对比图；
[0106] 图6为本发明的实施例提供的洱海子流域空间拓扑关系图；
[0107] 图7为本发明的实施例提供的弥苴河子流域空间拓扑关系矩阵；
[0108] 图8为本发明的实施例提供的洱海各子流域入湖污染贡献量空间分布图，(a)总氮,(b) 总磷；
[0109] 图9为本发明的实施例提供的洱海流域年降水量小波系数
实部系数等值线图；
[0110] 图10为本发明的实施例提供的洱海湖泊水质达标保证率P与流域污染负荷累积比例β 的对应关系(a.TN,b.TP)；
[0111] 图11为本发明的实施例提供的洱海流域总氮、总磷多级靶向治理目标识别结果；
[0112] 图12为本发明的目标水质—超标通量—入湖污染负荷映射关系结构图。
具体实施方式
[0113] 为了使本发明实施例的目的、技术方案和优点更加清楚，下面结合本发明提出的方法流 程图(图1)以及本发明实施例中的附图(图2～图11)，对本发明实施例中的技术方案进行 详细的说明。显然，所描述的实施例仅为本发明一部分实施例，而不是全部的实施例。
[0114] 本发明提供一种基于链式反应的湖泊型流域非点源污染多级靶向治理目标确定方法，该 方法中，克服了传统方法未考虑湖泊型流域大水体水质响应过程，忽略了子流域之间的水力 联系，以及靶向治理目标空间尺度识别的粗放化。在流域非点源污染靶向治理目标识别过程 中将流域水环境看成一个完整的系统来对待，充分考虑了污染物质从产生释放、输出入河、 迁移转化、再到水体水质响应这一完整的链式反应过程。基于流域分布式水文、非点源污染 模型SWAT依次得到子流域空间关系拓扑矩阵R、子流域污染输出负荷矩阵E、河网水系输 移影响关系矩阵F，从而进一步得到子流域入湖污染贡献量计算矩阵D；基于湖泊水体一维 完全混合降解模型，计算湖泊水体水质污染物响应系数ki；基于水文周期小波方法，计算湖 泊水环境评估周期T；基于灰色系统理论，构建目标水质—超标通量—入湖污染负荷映射关 系结构；基于上述计算结果，最终进行流域非点源污染多级靶向治理目标的识别工作(图1)。
[0115] 一种基于链式反应的湖泊型流域非点源污染多级靶向治理目标确定方法，首先进行数据 采集：收集流域气象、下垫面情况以及水文水质监测数据；然后进行多级靶向治理目标确定， 确定过程具体包括以下步骤：
[0116] 1)构建流域分布式水文、非点源污染模型：模型输入数据中包括流域地理高程、土地利 用类型、土壤类型以及坡度变化，模型以水文响应单元为基本计算网格，以子流域为统计单 元，子流域由多个水文响应单元组成，每个水文响应单元内具有相似的土地利用类型、土壤 类型和坡度；采用重要控制断面、子流域出口以及流域总出口的水文、水质监测数据进行模 型率定工作，获取适合于目标流域水文、营养物循环的关键参数；
[0117] 其中，选取流域分布式水文模型SWAT作为模拟流域水文、污染物循环的基础工具。
[0118] SWAT模型水量平衡方程如式(1)所示：
[0119]
[0120] 式中：SWt为土壤最终含水量，mm；SW0为土壤前期含水量，mm；t为时间步长，d；Rday为第i天降雨量；Qsurf为第i天的地表径流，mm；Ea为第i天的蒸发量，mm；Wseep为第i 天存在于土壤坡面底层的渗透量和侧流量，mm；Qgw为第i天地下水出流量，mm；
[0121] 泥沙输移方程如式(2)所示：
[0122] sed＝11.8(Qsurf·qpeak·areahru)0.56·Kusle·Cusle·Pusle·LSusle·CFRG (2)[0123] 式中：sed为土壤侵蚀量，t；Qsurf为地表径流量，mm/hm2；qpeak为峰值流量，m3/s；areahru为水文响应单元的面积，hm2；Kusle为土壤侵蚀因子；Cusle为植被覆盖和管理因子；
Pusle为保 持措施因子；LSusle为地形因子；CFRG为粗糙碎屑因子；
[0124] 氨氮硝化如式(3)所示：
[0125]
[0126] 式中：Nnit,ly为第ly层土层中 转化为 的氮含量，kg N/ha；frnit,ly为硝化过程中产生 的氮损失分数；frvol,ly为挥发过程中产生的氮损失分数；Nnit|vol,ly为第ly层土层硝化过程和挥发 过程共同引起的氮损失量，kg N/ha；
[0127] 硝酸盐反硝化如式(4)所示：
[0128]
[0129] 式中：Ndeniy,ly为反硝化过程所产生的氮损失量，kg N/ha；NO3ly为第ly层土层中的硝酸盐 含量，kg N/ha；βdenit为反硝化速率；γtmp,ly为第ly层土层中营养物循环温度因子；γsw,ly为 第ly层土层中营养物循环水分因子；orgCly为土壤层中有机碳含量，％；γsw,thr为营养物循环 中发生反硝化过程的水分因子阈值；
[0130] 残留物有机磷的矿化如式(5)所示：
[0131] Pminf,ly＝0.8·δntr,ly·orgPfrsh,ly (5)
[0132] 式中：Pminf,ly为残留物有机磷经过矿化作用转化的无机磷，kg N/ha；δntr,ly为残留物降解 系数；orgPfrsh,ly为第ly层土层中的残留物有机磷含量，kg/ha；
[0133] 溶解态无机磷与活性态有机磷之间的磷平衡方程如式(6)所示：
[0134]
[0135] 式中：Psol|act,ly为溶解态无机磷与活性态无机磷之间转化的磷含量，kg P/ha；Psolution,ly为溶 解态无机磷含量，kg P/ha；Pact,ly为活性态无机磷含量，kg P/ha；pai为磷可利用指数；当Psol|act,ly大于0时，磷从溶解态无机磷向活性态无机磷转化；当Psol|act,ly小于0时，磷从活性态无机磷 向溶解态无机磷转化。
[0136] 2)构建子流域空间拓扑关系矩阵R：相邻的子流域之间具有明确的上下游隶属关系；水 量和污染物质以子流域作为传输媒介；对于相邻的两个子流域Wi和Wj，Wi为紧邻Wj的上游 流域，其水量和污染物质由Wi转移到Wj的条件概率为1，反之，条件概率为0。
[0137] 对于n个子流域W1，W2，…，Wi，Wj，…，Wn，其空间拓扑矩阵可用n×n的矩阵进行 表示，Wi为紧邻Wj的上游流域，其水量和污染物质由Wi转移到Wj的条件概率为1，反之， 条件概率为0，其子流域空间拓扑关系矩阵数学表达式为：
[0138]
[0139]
[0140] 矩阵Rn×n的行向量之和为0或1；当行向量为1时，说明子流域有唯一确定的出口；当 行向量为0时，说明该子流域位于流域总出口处；在经过最多n次转移后，某一子流域Wi的水量和污染物质将从流域总出口流出，即：
[0141]
[0142] 式中：m为某一子流域Wi从源头到流域总出口所经历的最多子流域个数。
[0143] 3)构建子流域污染输出负荷矩阵E：根据已建立的流域分布式水文、非点源污染模型， 模拟计算子流域内的污染物质的产生、迁移和转化情况，某一子流域Wi的污染物输出量Ei是该子流域范围内所有水文响应单元污染物输出量的总和。
[0144] 子流域Wi的污染物输出量Ei是该子流域范围内所有水文响应单元污染物输出量的总和， 其计算公式如下：
[0145]
[0146] 式中：Ei为子流域Wi的污染物输出量；AHRU(i,j)为子流域Wi内HRUj的面积；CHRU(i,j)为 子流域Wi内HRUj的污染负荷输出强度；HRUj为水文响应单元；
[0147] 根据上述获取的不同子流域污染物输出量Ei，定义并建立各子流域的污染物输出量矩阵， 其数学表达式为：
[0148] E＝(E1,E2,···,Ei,···,En)T (15)
[0149] 4)构建河网水系输移影响关系矩阵F：根据各子流域的污染负荷输入及输出结果，采用 污染物通量相对变化值来表征某一段河网水系对于污染物输移的影响，河网水系输移影响系 数的取值范围在0到1之间，影响系数的数值越小，说明河网水系对于污染物输移影响的程 度越大，反之，河网水系对于污染物输移影响的程度越小。
[0150] 河道影响系数矩阵的数学表达式如下：
[0151]
[0152] 式中：Fi为子流域i的河道影响系数，其取值范围在0到1之间，河道影响系数表示河道 对污染物的削减或者贡献作用，本质上是污染物在河道中迁移转化的效率：
[0153]
[0154]
[0155]
[0156] 式中：Lout为子流域Wi出口断面的污染负荷量；Lin为子流域Wi入口断面的污染物负荷 量，其值为上游子流域Wi-1污染负荷的输出量Lout,Wi-1和子流域Wi的污染负荷输入量LWi之和。
[0157] 5)计算湖泊水体水质污染物响应系数k：采用水体一维完全混合降解模型获取湖泊水体 水质对于污染物质的响应系数k。
[0158] 步骤5)中，对于一个具有m个污染源的湖泊型流域，湖泊水体一维完全混合降解模型 的数学表达式如式(20)，所述污染源即入流湖泊：
[0159]
[0160] 式中：c为湖泊污染物浓度，mg/L；t为研究时段，d；V为湖泊水量，m3；qi为第i个污 染源的入湖流量，m3/d；ci为第i个污染源入湖污染物浓度，mg/L；Q为出湖水量，m3/d，kj为某一污染物的综合降解系数，1/d；
[0161] 对于一个长期存在且相对稳定的湖泊来说，若已知入湖水量、水质，出湖水量，湖泊水 量、水质的基础上，即可以求得某一污染物的综合降解系数，其数学表达式为：
[0162]
[0163] 6)构建子流域入湖污染贡献量矩阵D：基于上述建立的子流域空间拓扑关系矩阵R、子 流域污染输出负荷矩阵E、河网水系输移影响关系矩阵F以及湖泊水体水质污染物响应系数 k，计算各子流域入湖污染贡献量Di。
[0164] 子流域入湖污染贡献量计算数学表达式为：
[0165] Di＝(Rn×n·Fn×n)i·Ei·kj (22)
[0166] 式中：Di为某一子流域Wi的入湖污染物贡献量kg，i≤n；Rn×n为子流域空间拓扑关系矩 阵；Fn×n为河道影响系数矩阵；Ei为某一子流域污染输出量，kg；kj为某一污染物的综合降 解系数，1/d。
[0167] 7)确定湖泊水环境评估周期T：采用水文小波分析方法确定湖泊水环境评估周期T。
[0168] 引入基本小波函数Ψ(t)，通过小波方差图，可以确定一个水文序列中存在的主要时 间尺度，即主周期。
[0169] 8)构建湖泊水环境状态变化、湖区水体污染物超标通量变化和流域污染负荷变化一一映 射关系结构：分别通过湖泊目标水质达标保证率P、湖泊污染物超标积累通量α及流域污染 负荷累积比例β表示。
[0170] 目标水质达标保证率P的数学表达式为：
[0171]
[0172] 式中：T达标为湖区水质达到目标水质的时间，T评估周期为评估周期时长；
[0173] 湖泊污染物超标累积通量α的数学表达式为：
[0174]
[0175] 式中：C观测为湖区水质观测值，mg/L；C目标为水质目标值，mg/L；为了保险起见，以 年内月份最差的水质浓度作为该年湖泊污染物超标量的计算依据；V代表水位为典型代表水文年 下特征水位所对应的湖泊水量，m3；α的单位为kg。
[0176] 与之相对应的，流域污染负荷累积比例β的数学表达式为：
[0177]
[0178] 式中：α表示湖泊污染物超标累积通量，kg；kj为某一污染物的综合降解系数，1/d；M为流 域非点源治理措施实施效率系数；L总负荷量为典型代表水文年的流域入湖总污染负荷量，kg。。 β单位为％。
[0179] 9)识别不同水质达标率下的流域非点源污染多级靶向治理目标：根据步骤6)得到各 子流域污染贡献量Di，为了消除子流域面积对于计算结果的影响，采用污染贡献量强度Ii对 各子流域进行排序，并绘制子流域污染贡献量累积曲线CD和子流域面积累积曲线CA；采用 步骤8)得到的不同湖泊水质达标率下所对应的流域污染负荷累积比例作为流域非点源污染 靶向治理目标识别的划分依据，获取为达到某一湖泊水质达标率时需要优先控制的子流域位 置、面积信息。
[0180] 实施例1
[0181] 本发明实施例将一种基于链式反应的湖泊型流域非点源污染多级靶向治理目标确定方法 应用在
云南省洱海流域(图2)，采用如下步骤进行实施：
[0182] (1)洱海流域分布式水文、非点源污染模型的构建。基础
数据库包括空间数据库(地理 高程、水下地形、土地利用、土壤类型)、属性数据库(土地利用类型属性、土壤类型物化属 性、气象要素数据、水文水质数据)以及其他数据库(作物管理措施、社会经济数据、流域 管理措施)。在基础数据库的基础上，选取典型入湖河流及重要控制断面的水文、水质监测资 料进行模型率定，获取适合洱海流域的水文、营养物循环关键参数，其模拟结果与实测数据 的对比图如图3所示。
[0183] (2)子流域空间拓扑关系矩阵R。根据子流域上下游的水力联系，建立子流域空间拓扑 关系，图4展示了洱海流域各条入湖河流子流域的上下游空间拓扑关系图。基于流域分布式 水文、非点源污染SWAT模型，洱海流域共划分为87个子流域，73号和74号子流域位于洱 海的出海口，对洱海湖泊的水环境影响很小，未纳入本次研究范围。其中，北部洱源平地河 流区涉及到弥苴河、罗时江和永安江三条主要入湖河流，共划分有27个子流域；西部大理高 山河流区涉及到苍山十八溪、棕树河和21条主要入湖沟渠，共划分有40个子流域；南部凤 仪中山河流区涉及到波罗江和白塔河两条主要入湖河流，共划分有15个子流域；东部海东低 山河流区涉及到凤尾箐、玉龙沟、南村河和下和箐四条主要入湖河流，共划分有4
2
个子流域。 子流域平均面积为26.82km ，占洱海总流域面积的1.15％，足以能反映洱海流域复杂的气候 和下垫面情况。考虑到洱海子流域众多，需要建立87×87阶的矩阵进行表示，受限于篇幅， 以洱海流域最大的入湖河流弥苴河作为示例来展示子流域水力联系关系转化为空间拓扑关系 矩阵，其结果如图4所示。
[0184] (3)子流域污染输出污染负荷矩阵E。对SWAT模型的输出文件SWATOutput.mdb进行 汇总(图5)，提取rch数据表中的TOT_Nkg和TOT_Pkg作为各子流域的污染物输出数据。
[0185] (4)河网水系输移影响关系矩阵F。对SWAT模型的输出文件SWATOutput.mdb进行汇 总(图5)。提取rch数据表中的子流域进口有机氮(ORGN_INkg)、有机磷(ORGP_INkg)、 硝酸盐(NO3_INkg)、氨氮(NH4_INkg)、亚硝酸盐(NO2_INkg)、无机磷(MINP_INkg) 的作为式(17)中Lin的输入数据。提取rch数据表中的子流域出口有机氮(ORGN_OUTkg)、 有机磷ORGP_OUTkg)、硝酸盐(NO3_OUTkg)、氨氮(NH4_OUTkg)、亚硝酸盐 (NO2_OUTkg)、无机磷(MINP_OUTkg)的作为式(17)中Lout的输入数据。根据式(17) 计算得到各子流域的河网水系输移影响关系矩阵。
[0186] (5)湖泊水体水质污染物响应系数k。根据洱海流域不同年份的降水量，选取典型的丰 水年、平水年和枯水年。根据典型代表年的洱海湖泊入湖水量水质、出湖水量水质以及湖区 水质监测站的水质监测数据，基于式(21)，计算得到洱海湖泊水体对于总氮、总磷的水质响 应系数。
[0187] (6)子流域入湖污染贡献量矩阵D。以子流域为统计单元，基于上述建立的子流域空间 拓扑关系矩阵R、子流域污染输出污染负荷矩阵E、河网水系输移影响关系矩阵F以及湖泊 水体水质污染物响应系数k，计算和统计洱海各子流域多年平均总氮、总磷入湖污染贡献量， 其空间分布见图6。每个子流域的污染物质均经历了源强产生、子流域输出、河网水系的迁 移转化以及水体水质响应的链式反应过程。总体上来看，总氮、总磷入湖污染贡献量的空间 分布较为相似，总氮入湖污染贡献量介于0kg～70490kg之间，总磷入湖污染贡献量在 0kg～4554kg之间。整体上表现为北部洱源平地河流区、西部大理高山河流区的入湖污染贡献 量大于南部凤仪中山河流区、东部海东低山河流区的入湖污染贡献量。
[0188] (7)湖泊水环境评估周期T。采用水文周期小波分析方法对洱海流域1951年～2017年近 67年的年降水量进行周期分析。可以看出，基于小波变换识别出1951年～2017年洱海流域年 降水量具备三个主周期，分别为55a、10a和35a(图7)。考虑到湖区水质监测资料的完整性， 选取10a作为湖泊水环境状态的评估周期，即2008年～2017年。
[0189] (8)构建目标水质—超标通量—入湖负荷映射关系结构。根据《水功能区划分标准》 (GB/T50594-2010)，洱海湖泊被划分为苍山洱海自然保护区，其水功能区水质要求为Ⅱ类水， 对应的总氮、总磷浓度标准分别为0.5mg/L和0.025mg/L。从全湖平均水质状况来看，洱海多 年全湖平均总氮浓度介于0.53mg/L～0.83mg/L之间，全湖平均总磷浓度介于 0.0285mg/L～0.0433mg/L之间。但总体上来看，从最严格角度来看，近10年的洱海全湖平均 水质均未达到Ⅱ类水质要求。为了便于相关部门进行湖区水环境治理，本文将湖区水质达标 率统一划分为三类类，即湖区水质达标保证率达到50％、湖区水质达标保证率达到80％和湖 区水质达标保证率达到100％。考虑到洱海流域非点源污染治理削减效果，本文将洱海流域非 点源治理措施实施效率系数γ设定为0.85。洱海湖泊2008年～2017年水质达标保证率P和流 域污染负荷累积比例β的对应关系如图8所示。
[0190] (9)流域非点源污染多级靶向治理目标识别。以典型平水年2017年为例，根据各子流 域入湖污染贡献强度对各子流域进行排序，图9和表1-2展示了针对全湖平均水质状态的洱 海流域总氮、总磷多级靶向治理目标识别结果。在总氮方面，入湖污染贡献负荷强度在 0.04kg/ha～66.44kg/ha之间。其中，有17个子流域为一级靶向治理目标，其累积入湖污染负 荷和累积面积百分比分别为48.03％和6.14％；有8个子流域为二级靶向治理目标，其累积入 湖污染负荷和累积面积百分比分别为57.75％和9.33％；有11个子流域为三级靶向治理目标， 其累积入湖污染负荷和累积面积百分比分别为79.86％和21.49％；其余有45个子流域为四级 靶向治理目标。在总磷方面，入湖污染贡献负荷强度在0.001kg/ha～3.900kg/ha之间。其中， 有21个子流域为一级靶向治理目标，其累积入湖污染负荷和累积面积百分比分别为51.93％ 和8.74％；有16个子流域为二级靶向治理目标，其累积入湖污染负荷和累积面积百分比分别 为75.60％和22.99％；有14个子流域为三级靶向治理目标，其累积入湖污染负荷和累积面积 百分比分别为87.39％和34.16％；其余有30个子流域为四级靶向治理目标。
[0191] 总体上可以看出，在同等的累积子流域面积上，其累积入湖污染量也不尽相同，80％的 入湖氮磷污染负荷量集中在流域30％以内的区域。值得注意的是，随着湖泊水质达标保证率 的升高，子流域累积入湖污染负荷量的
曲率逐渐变小，而子流域面积累积曲线的曲线逐渐增 大，说明实现的湖泊水质达标保证率越高，其需要治理的流域面积会显著增大，其难度也会 相应增加。具体的，在总氮、总磷的累积入湖污染负荷量分别小于60％和
55％时，其子流域 累积面积曲线的曲率处于相对较小的变化范围。因此，对于总氮来说，将湖泊水质达标保证 率从10％提升至80％时较为经济(控制9.33％的流域面积)，对于总磷来说，将湖泊水质达标 保证率从10％提升至50％时较为经济(控制8.74％的流域面积)。
[0192] 附表1为本发明的实施例提供的洱海各子流域总氮污染多级靶向治理目标识别结果(平 水年)；
[0193] 附表2为本发明的实施例提供的洱海各子流域总磷污染多级靶向治理目标识别结果(平 水年)。
[0194] 上述的实施例仅是本发明的部分体现，并不能涵盖本发明的全部，在上述实施例以及附 图的基础上，本领域技术人员在不付出创造性劳动的前提下可获得更多的实施方式，因此这 些不付出创造性劳动的前提下获得的实施方式均应包含在本发明的保护范围内。
[0195] 附表1洱海各子流域总氮污染多级靶向治理目标识别结果(平水年)
[0196]
[0197]
[0198]
[0199] 附表2洱海各子流域总磷污染多级靶向治理目标识别结果(平水年)
[0200]
[0201]
[0202]
