基于水质-水位二元响应关系的适宜生态水位预测方法

1.本发明涉及基于水质-水位二元响应关系的适宜生态水位预测方法，属于水位预测技术领域。


背景技术：

2.随着全球变暖和社会经济的快速发展，国内众多湖泊出现了面积萎缩、水质恶化、生态退化等问题。而水位是反映湖泊水文情势和湖泊规划管理的重要指标，其时空变化过程会对湖泊生态系统产生深远的影响，是解决上述问题的关键所在。关于湖泊生态水位的研究分为两个方面：最小生态水位和适宜生态水位。最小生态水位是满足湖泊生态系统基本需求的水位红线，其往往可在短时期内保护湖泊生态系统功能和结构不受到严重破坏，但长期地维持会导致湖泊流域生态系统退化甚至崩溃；并且最小生态水位多为一固定的水文特征值，难以反映湖泊系统在不同时间节点对水位的要求。因此，确定维持湖泊生态系统可持续健康发展的适宜生态水位成为了湖泊生态研究的新热点。
3.目前，国内外学者多借鉴iha-rva法计算湖泊天然水位情势，以此确定湖泊适宜生态水位需求。但iha-rva法侧重于还原历史资料中的天然水位情势，展现人类干预较少情况下湖泊生态系统对于水量的需求。大型过水型湖泊往往承担着防洪、供水和航运等多种功能，剧烈的人类活动改变了国内大型过水型湖泊的生态环境，污染负荷的增加使水质状况出现了不同程度的恶化。重金属富集和蓝藻爆发等污染现象时有发生，破坏了湖泊生态系统健康的可持续发展。仅依靠iha-rva方法计算的适宜生态水位难以全面地反映生物群落健康发展对生存环境的需求，有必要将水质因素纳入过水型湖泊适宜生态水位推求的考虑范围。
4.目前，国内学者针对过水型湖泊水位-水质相关性开展了大量研究，但考虑过水型湖泊水质需求的适宜生态水位研究还较为少见，管理部门多使用最小生态水位作为过水型湖泊的生态需水指标。同时过水型湖泊出入湖水量大的特点，也使现有水质调控的一些成熟经验如特定时期增加入湖水量改善湖区水质等的实践应用效果不佳。


技术实现要素：

5.本发明所要解决的技术问题是克服现有技术的缺陷，提供基于水质-水位二元响应关系的适宜生态水位预测方法。
6.为达到上述目的，本发明提供基于水质-水位二元响应关系的适宜生态水位预测方法，包括：(1)根据水文要素和水质要素的变化趋势，划分年内时段；
7.(2)根据滞留污染物总量、纳污能力和水位-水质经验公式，确定满足设定水质目标的水位区间和逐月水位变动率；
8.(3)采用iha-rva法计算天然水位情势；
9.根据水位区间和水位-水质响应关系修正天然水位情势，得到满足设定水质目标的适宜生态水位。
10.优先地，步骤(1)，具体包括：
11.获取一段时间内湖泊的水文要素和水质要素，计算获得多年月平均出湖水量、年平均出湖水量、多年月平均入湖污染物浓度和多年月平均入湖污染物总量；
12.根据水文要素和水质要素的年内变化趋势，划分年内时段。
13.优先地，水文要素包括湖泊水位、入湖流量和出湖流量；
14.水质要素包括湖泊污染物浓度和入湖污染物浓度；
15.年内时段包括水位波动时段、水位上升时段和水位下降时段。
16.优先地，步骤(2)，包括：
17.选取发生频率为x和y的多年月水位变化率作为逐月水位变化率的阈值，x大于y，x和y为常数，采用湖泊水位变动率衡量滞留污染物总量：
18.δhz
′i＝p
x
(δz
′
)
ꢀꢀ
(1)，
19.δlz
′i＝py(δz
′
)
ꢀꢀ
(2)，
20.式中，δhz
′i为逐月水位变化率上限，δlz
′i为逐月水位变化率下限；p
x
(δz
′
)为发生频率x的湖泊多年月水位变化率，py(δz
′
)为发生频率y的湖泊多年月水位变化率。
21.优先地，步骤(2)，确定满足设定水质目标的水位区间，包括：
22.采用dillon模型计算湖泊的纳污能力；
23.基于多年月平均出湖水量和多年月平均入湖污染物总量，计算逐月适宜水位下限：
[0024][0025][0026]
qa＝q
×
12/n
ꢀꢀ
(5)
[0027]
式中，w为湖泊中总氮的纳污能力；ρs为总氮控制目标浓度；z为湖泊平均深度；qa为湖泊的年平均出湖水量；a为湖泊面积；v为湖泊库容；r
p
为总氮在湖泊中的滞留系数；q为计算月的出湖水量；n为计算月的月数；
[0028]
基于选用的一段时间的湖泊水位和湖泊污染物浓度，绘制不同时段内的湖泊水位-水质散点图，拟合水位-水质经验公式；
[0029]
基于水位-水质经验公式，计算获得满足设定水质目标的水位区间。
[0030]
优先地，步骤(2)，确定满足设定水质目标的逐月水位变动率，包括：
[0031]
采用月平均水位z
ave，i
和若干日的平均日水位变化率δzj作为适宜生态水位的计算指标，i＝1～12，j＝1～365；
[0032]
选取发生频率为x和y的计算指标作为rva阈值，确定不同时期湖泊对满足设定水质目标的逐月水位变动率。
[0033]
优先地，步骤(3)，采用iha-rva法计算天然水位情势；
[0034]
根据水位区间和水位-水质响应关系修正天然水位情势，得到适宜生态水位，包括：
[0035]
水位波动时段中，利用满足设定水质目标的水位区间替换设定的不满足设定水质目标的天然水位情势；
[0036]
基于水位-水质响应关系，对水位波动时段的天然水位情势进行逐月水位阈值修正；
[0037]
水位-水质响应关系的计算公式为：
[0038][0039][0040]
式中：hz
′i为修正后水位波动时段的第i月适宜生态水位上限；lz
′i为修正后水位波动时段的第i月适宜生态水位下限；hz
wq
为满足设定水质目标的水位上限；lz
wq
为满足设定水质目标的水位下限；hzi为修正前水位波动时段的第i月适宜生态水位上限；lzi为修正前水位波动时段的第i月适宜生态水位下限。
[0041]
优先地，水位上升时段和水位下降时段中，利用设定水质目标的水位上限和设定水质目标的水位下限限制水位总变动幅度，基于逐月水位变化率上限、逐月水位变化率下限和计算获得的湖泊的纳污能力，分配水位变动过程；
[0042]
水位变动过程分配的计算公式如下：
[0043][0044][0045]
式中：hz
′k为修正后第k月适宜生态水位上限；lz
′k为修正后第k月适宜生态水位下限；hzk为修正前第k月适宜生态水位上限；lzk为修正前第k月适宜生态水位下限；hz0为水位变动前一月修正后适宜生态水位上限；lz0为水位变动前一月修正后适宜生态水位下限；lz
p
为基于水位-水质经验关系的计算公式计算出的修正后水位波动时段的第p月适宜生态水位上限；lzr为基于水位-水质经验关系的计算公式计算出的修正后水位波动时段的第r月适宜生态水位下限；lz
p
为基于纳污能力计算出的水位下限；j为水位变化期间的总月数；k＝1，2，
…
，j-1。一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现上述任一项所述方法的步骤。
[0046]
优先地，一种计算机可读存储介质，其上存储有计算机程序，该计算机程序被处理器执行时实现上述任一项所述方法的步骤。
[0047]
本发明所达到的有益效果：
[0048]
1、本发明从滞留污染物总量、纳污能力和水位-水质经验公式三个方面，分别量化不同时期内的水位-水质响应关系，得到了满足湖泊水质目标的水位区间。
[0049]
2、本发明在适宜生态水位设计过程中引入水位-水质响应关系，增强了设计的适宜生态水位的适用性，满足了湖泊生态系统对于水质的需求。
[0050]
3、本发明为大型过水型湖泊的生态调度和水资源管理提供参考和依据，具有很强的实用性和广泛的适用性。
附图说明
[0051]
图1为本发明实施例二的流域图；
[0052]
图2为本发明修正水位波动时段的适宜生态水位的示意图；
[0053]
图3为本发明修正水位上升时段的适宜生态水位的示意图；
[0054]
图4为本发明修正水位下降时段的适宜生态水位的示意图；
[0055]
图5为洪泽湖入湖水量、出湖水量和湖泊水位的逐月变化趋势图；
[0056]
图6为洪泽湖五日平均水位变动率的示意图；；
[0057]
图7为洪泽湖入湖水量和水位的关系图；
[0058]
图8为水质因素逐月变化趋势图；
[0059]
图9为平水期水位和水质站点达标比例的散点图；
[0060]
图10为泄水期水位和水质站点达标比例的散点图；
[0061]
图11为蓄水前期水位和水质站点达标比例的散点图；
[0062]
图12为蓄水后期水位和水质站点达标比例的散点图；
[0063]
图13为本发明iha-rva法计算的适宜生态水位的示意图；
[0064]
图14为本发明iha-rva法计算的适宜逐日水位变动率图；
[0065]
图15为本发明修正后适宜生态水位的区间图。
具体实施方式
[0066]
以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
[0067]
实施例一
[0068]
基于水质-水位二元响应关系面向过水型湖泊推求适宜生态水位的方法，包含步骤：
[0069]
(1)根据水文要素和水质要素的变化趋势，划分年内时段；
[0070]
(2)从滞留污染物总量、纳污能力和水位-水质经验公式三个方面，分别量化不同时期内的水位-水质响应关系，并基于此确定满足设定水质目标的水位区间a和逐月水位变动率；
[0071]
(3)采用iha-rva法计算湖泊的天然水位情势，并根据水位区间a修正天然水位情势，得到满足设定水质目标(如v类水质)的适宜生态水位。
[0072]
所述步骤(1)根据水文要素和水质要素的变化趋势划分年内时段：
[0073]
采用湖泊30年以上的水文要素和5年以上的水质要素获得多年月平均出湖水量、多年月平均入湖污染物浓度和多年月平均入湖污染物总量，分别分析水文要素和水质要素的年内变动趋势，并基于此划分年内时段；
[0074]
水文要素包括湖泊水位、入湖流量和出湖流量；
[0075]
水质要素包括湖泊污染物浓度和入湖污染物浓度；
[0076]
年内时段包括水位波动时段、水位上升时段和水位下降时段。
[0077]
所述步骤(2)，从滞留污染物总量、纳污能力和水位-水质经验公式三个方面，分别量化不同时期内的水位-水质响应关系，并基于此确定满足水质目标的水位区间a，包括：
[0078]
(21)采用湖泊水位变动率衡量滞留污染物总量，选取预设发生频率为75％和25％
的多年月水位变化率作为逐月水位变化率的阈值，计算公式见下：
[0079]
δhz
′i＝p
75％
(δz
′
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0080]
δlz
′i＝p
25％
(δz
′
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0081]
式中：δhz
′i和δlz
′i分别为湖泊逐月水位变化率上限和下限，m/月；p
75％
(δz
′
)和p
25％
(δz
′
)为75％和25％的发生频率的湖泊多年月水位变化率，m/月。
[0082]
(22)采用dillon模型计算湖泊的纳污能力，并基于多年月平均出湖水量和多年月平均入湖污染物总量，计算逐月适宜水位下限，计算公式为：
[0083][0084][0085]
qa＝q
×
12/n
ꢀꢀ
(5)
[0086]
式中：w为湖泊中总氮的纳污能力，t/a；ρs为总氮控制目标浓度，mg/l；z为湖泊平均深度，m；qa为湖泊的平均年出湖水量，m3；a为湖泊面积，m2；v为湖泊库容，m3；r
p
为总氮在湖泊中的滞留系数，a-1
；q为计算月的出湖水量，m3；n为计算月的月数。
[0087]
(23)采用5年以上的湖泊水位和污染物浓度，绘制不同时段内的湖泊水位-水质散点图，并拟合水位-水质经验公式；
[0088]
基于水位-水质经验公式，计算获得满足设定水质目标的水位区间。
[0089]
所述步骤(3)，采用iha-rva法计算湖泊的天然水位情势，并根据水位区间a修正天然水位情势，得到满足设定水质目标的适宜生态水位，具体包括：
[0090]
(1)基于洪泽湖多年逐日水位数据，采用月平均水位(z
ave，i
，i＝1～12，m)和五日平均日水位变化率(δzj，j＝1～365，m)作为适宜生态水位计算指标，选取发生频率为25％和75％的计算指标作为rva阈值，初步确定不同时期湖泊生态系统对湖泊水体大小与水位变化速率的需求。
[0091]
(2)根据湖泊水位变动趋势将划分的年内时段分为水位波动时段、水位上升时段和水位下降时段。如图2所示，水位波动时段中以计算获得满足设定水质目标的水位区间直接替换不满足设定水质目标的适宜生态水位阈值，基于水位-水质响应关系对水位波动时段的iha-rva计算结果进行逐月水位阈值修正；
[0092]
水位-水质响应关系的计算公式如下：
[0093][0094][0095]
式中：hz
′i为修正后水位波动时段的第i月适宜生态水位上限，m；lz
′i为修正后水位波动时段的第i月适宜生态水位下限，m；hz
wq
为满足设定水质目标的水位上限，m；lz
wq
为满足设定水质目标的水位下限，m；hzi为修正前水位波动时段的第i月适宜生态水位上限，m；lzi为修正前水位波动时段的第i月适宜生态水位下限，m。
[0096]
对于水位上升时段和水位下降时段中，应采用式(6)-(7)，通过满足设定水质目标的水位上限和满足设定水质目标的水位下限限制水位总变动幅度，再按五日平均的逐月水
位变化率上限、逐月水位变化率下限和计算获得的湖泊纳污能力，分配水位变动过程，如图3和4所示。水位变动过程分配的计算公式如下：
[0097][0098][0099]
式中：hz
′k为修正后第k月适宜生态水位上限，m；lz
′k为修正后第k月适宜生态水位下限，m；hzk为修正前第k月适宜生态水位上限，m；lzk为修正前第k月适宜生态水位下限，m；hz0为水位变动前一月修正后适宜生态水位上限，m；lz0为水位变动前一月修正后适宜生态水位下限，m；lz
p
为基于水位-水质经验关系的计算公式计算出的修正后水位波动时段的第p月适宜生态水位上限，m；lzr为基于水位-水质经验关系的计算公式计算出的修正后水位波动时段的第r月适宜生态水位下限，m；lz
p
为基于纳污能力计算出的水位下限，m；j为水位变化期间的总月数；k＝1，2，
…
，j-1。
[0100]
采用iha-rva法计算天然水位情势，采用dillon模型计算湖泊的纳污能力，依据湖泊水位和湖泊污染物浓度的历史数据拟合水位-水质经验公式，均属于现有技术，本实施例不再详细阐述。
[0101]
实施例二
[0102]
本发明选取洪泽湖作为实施案例，洪泽湖地处淮河中下游结合部，汇水面积达15.8km2，是中国第四大淡水湖。作为过水型湖泊，洪泽湖吞吐性强，换水周期短。自上世纪90年代以来，洪泽湖污染逐渐严重，水质恶化，近年来虽有好转，但仍处于轻度富营养化状态。目前，洪泽湖水质多为劣v类，总氮是主要污染物，不考虑总氮的情况下洪泽湖水质可达iii类标准。洪泽湖水体主要依靠地表径流补给，其主要入湖河流为淮河、新汴河、怀洪新河、池河、新濉河、老濉河和徐洪河，主要出湖河流为入江水道、苏北灌溉总渠和入海水道。其中淮河是洪泽湖水量的主要来源，也是入湖总氮的主要来源。洪泽湖地理位置及出入湖河流和水质监测点分布如图1所示。
[0103]
(1)根据水文要素和水质要素的变化趋势划分年内时段
[0104]
考虑到洪泽湖的自身特点，本发明选取总氮为洪泽湖污染物，并以v类水质为水质目标进行计算。采用1988-2018年的水文数据和2013-2018年的水质数据计算平均值并绘制变化趋势图，如图5所示。多年逐月入湖水量变化范围为8.0～86.5亿m3，全年可分为三个水量变化时期：1-6月洪泽湖入湖水量较少，总体呈波动趋势；7月份入湖水量达到峰值，随后8-9月呈下降趋势，但总量较大，7-9月占全年入湖水量的65.9％；10-12月入湖水量大幅下降，12月的入湖水量不足10月的一半。按五日水位变化率全年可以分为四个时期：1-4月洪泽湖水位在13～13.14m之间小幅波动；5-6月由于汛期的防汛需求，洪泽湖水位迅速下降；7月后随着汛期来水湖泊水位迅速上涨；汛期结束后，10-12月水位缓慢抬升。如图6所示为五日平均水位变动率，如图7所示，洪泽湖入湖水量和水位的关系为时间序列循环曲线，可被分为4个具有明显特征的阶段。
[0105]
水质因素变化趋势图如图8所示，根据入湖总氮浓度和入湖总氮总量，可将洪泽湖污染负荷年内变化分为四个时期。洪泽湖逐月平均总氮入湖浓度范围为1.96～3.75mg/l，年内变化为：1-4月浓度在3.37-3.75mg/l之间波动；5-6月浓度呈小幅下降趋势；7-9月由于汛期入湖水量迅速增加，造成总氮浓度大幅下降，9月仅为1.75mg/l；10-12月总氮浓度出现大幅反弹。洪泽湖多年月平均总氮入湖总量变化范围为3170.44～12301.92t，汛期入湖量远大于非汛期，年内变化趋势以七月为界主要呈现单峰变化，先增后减，10月出现小幅反弹。
[0106]
基于洪泽湖水文情势和总氮负荷年内变化趋势(如下表)，将全年分为四个时期：平水期(1至4月)、泄水期(5至6月)、蓄水前期(7至9月)和蓄水后期(10至12月)，分别进行水位-水质响应关系的量化。
[0107][0108]
(2)从滞留污染物总量、纳污能力和水位-水质经验公式三个方面分别量化不同时期内的水位-水质响应关系，并基于此确定满足水质目标的逐月水位变动率和水位区间a
[0109]
(1)计算的洪泽湖适宜逐月水位变化率变化范围为-0.79m～1.12m如下所示。
[0110][0111][0112]
(2)以v类水质标准，根据dillion模型计算得出的洪泽湖不同时段内纳污能力和水位下限如下所示，其中平水期11.08m，泄水期12.39m，蓄水后期12.9m。
[0113][0114]
(3)采用洪泽湖2013-2018年平水期水位、泄水期水位、蓄水前期水位、蓄水后期水位和水质站点达标比例绘制不同时期的散点图，如图9-图12所示。蓄水前期水质有随水位升高而变好的趋势，但相关性较差，而余下各时段的水质均与水位呈显著负相关关系。去除
部分水位变动趋势异常数据后，采用spss软件对这三个时期的水位-水质站点达标比例进行相关性计算和经验关系公式拟合，结果如下所示。
[0115][0116]
(4)综合纳污能力和水位-水质经验公式计算出的水位要求，得到满足设定水质目标的水位区间a为平水期低于12.99m，泄水期12.39～12.63m，蓄水后期12.90～13.04m。
[0117]
(3)基于iha-rva法计算湖泊天然水位情势，并根据水位区间a修正天然水位情势，得到满足设定水质目标的适宜生态水位
[0118]
(1)基于iha-rva计算得出的适宜生态水位阈值和适宜逐日水位变动率分别见图13和图14，洪泽湖适宜生态水位变化范围为11.97～13.41m，适宜逐日水位变化范围为-0.06m～0.05m。如图13所示，非汛期洪泽湖水位变化幅度较小，汛期洪泽湖水位呈现先降后升的趋势；如图14所示，洪泽湖的大幅度水位变动主要发生在5-8月，其中5-6月呈下降趋势而7-8月呈上升趋势，年内其余时段的水位波动幅度较小。
[0119]
(2)根据水位-水质经验公式和纳污能力要求修正的逐月适宜生态水位区间如图15所示，各月水位区间为：12.92～12.99m、12.79～12.99m、12.84～12.99m、12.86～12.99m、12.71～12.89m、12.39～12.63m、11.97～12.93m、12.50～13.07m、12.65～13.26m、12.90～13.04m、12.90～13.04m、12.90～13.04m。
[0120]
本技术是参照根据本技术实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0121]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0122]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0123]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。