滨海地区基于避咸蓄淡模式的供水安全临界流量研究——以粤港澳大湾区珠澳供水为例*

涂新军,吴海鸥,陈晓宏,林凯荣,赖荣标

(1:中山大学土木工程学院水资源与环境研究中心， 广州 510275)

(2:广东省华南地区水安全调控工程技术研究中心，广州 510275)

(3:南方海洋科学与工程广东省实验室，珠海 519000)

感潮河流淡水资源开发利用最大制约因素是咸潮入侵[1]. 咸潮入侵是指入海口高盐度水体上溯，造成河道水体含氯度超标现象，其受到河口潮流涨落作用、入海径流冲淡作用和风力扰动等多方面影响[2]. 河口三角洲咸潮入侵研究受到普遍关注[3]，如发生在我国的闽江口[4]、长江口[5-6]和珠江口[7-8]等. 随着当地社会经济快速发展，区域用水持续增长，咸潮影响下供水安全问题日益突显[9]，如珠江口珠海市咸潮影响区的水资源优化配置研究[10]和长江口上海市多个取水口供水安全问题[11]. 当枯水期上游来水量偏少时，感潮河段取水口含氯度在咸潮影响下部分时段超标，不能达到取水标准，从而威胁城市供水安全. 通过流域水量调度增加枯水期下泄流量，成为保障感潮河段取水安全的措施之一. 但枯水期可利用水资源有限，如何确定合理的上游来水压咸临界流量，成为咸潮影响下感潮河段供水安全关键问题之一.

在长江河口，有人提出了保障长江口供水安全的月控制临界流量含义和标准，并结合水动力模型分析了长江口咸潮入侵期的临界流量过程[1]；基于径流和潮汐的共同影响，提出了南渡河口临界水位线和淡水保证率概念[12]；运用水动力和咸水入侵模型，定量分析上游来水对咸度入侵边界线影响并确定压咸流量[5]；综合多种盐度预测经验模型，确定了长江口南支上段压咸临界流量[13]. 在珠江河口，有人基于径流-咸潮响应模式，分析磨刀门水道供水的压咸流量和压咸时机[14]；通过探讨磨刀门地区盐水入侵距离与潮差的响应关系，得到增大压咸流量的最佳时期[7]. 前人对于压咸临界流量的研究中，主要采用已有短期咸度观测资料，与上游来水流量进行统计分析，较少考虑当地水库供蓄能力对临界流量的影响. 在分析长江口临界流量时，有研究认为水库由于其多目标性对压咸调控有限[6]. 如何基于上游来水进行感潮河段咸度预测，定量构建上游来水与超标过程之间的联系，以及避咸取水如何与当地水库调度结合，是保障当地供水安全管理的迫切需求.

珠江三角洲是我国粤港澳大湾区核心区域，磨刀门水道是珠海及澳门的重要水源地，枯水期咸潮入侵已经威胁到了当地供水安全[15]. 自2005年以来，通过西江上游骨干水库群水量调度，加大了枯水期西江下游来水量，对于抑制珠江河口磨刀门水道的咸潮上溯效果显著[16-17]. 随着大藤峡水利枢纽工程投入运行及与骨干群水库的压咸调度，提出了当枯水期咸潮上溯影响严重时，梧州断面下泄流量调度目标为2100 m3/s，能有效地保障珠江三角洲河口区域供水安全[18]. 但是一方面，近十年来磨刀门水道主要取水口，短时间尺度(时均值或日均值)咸度超标问题依然存在. 另一方面，随着粤港澳大湾区发展规划的实施，区域经济高速增长及城市群发展模式对区域供水安全提出了更高要求[19]. 近些年来，当地水资源规划提出了“江水为主、库水为辅、江库联动、江水补库、库水调咸”的供水模式. 上游来水流量及对应的避咸取水时段，也成为了当地河库联合调度供水关键问题.

因此，本文研究拟构建避咸蓄淡供水模型，基于BP神经网络的咸潮预测模块，模拟不同来水过程下咸潮超标情况；通过河库联合供水调度模块，分析不同需水情景下长系列上游来水日过程的供需平衡状况；通过供水安全分析模块，识别上游来水的风险流量及临界流量，为依赖感潮河段为水源地的区域供水安全管理提供科学依据. 以面向粤港澳大湾区珠海及澳门的珠江三角洲磨刀门水道取供水为例，结合当地水资源规划，分析不同需水情景下当地河库联合供蓄的上游临界流量.

1 研究方法

1.1 避咸蓄淡供水模型

当感潮河流取水口的咸度超标成为常态，且又无法寻找到替代水源时，如何充分利用感潮河流适宜的取水时机(咸度小于水质标准限值的取水时间)和水库等蓄水工程，通过避咸蓄淡及河库联合供水，成为了沿海城市(如珠江三角洲的珠海市)保障供水安全的重要措施[20]. 其中，在城市供水过程中，取水口超标时间和水库可调蓄能力，成为当地供水安全保障的关键参数. 为此，避咸蓄淡供水模型(withdrawal avoiding saltwater model，简称WASM)，耦合了咸度预测模块、河库联合供水调度模块和供水安全分析模块，如图1.

图1 感潮河流避咸蓄淡供水模型

通过咸度预测模块，建立上游来水量与取水口超标时间的联系，并根据上游来水日过程，预测取水口日超标时间. 通过河库联合供水调度模块，基于特定需水量，分析某一枯水期上游来水下的取水、供水、蓄水及缺水等日过程. 通过供水安全分析模块，基于是否缺水以及设定阈值下水库剩余水量判断，识别上游来水是否属于风险流量，再结合设定阈值下总缺水量判断，从风险流量中识别供水安全的临界流量.

1.2 咸度预测和超标时间

感潮河流咸度变化过程受多重因素影响，主要包括河口潮流涨落作用海水入侵、上游径流冲淡作用和风力对河口纵向净环流的扰动作用等[21]. 因此，可建立咸度因子S与其主要影响因子(潮汐因子H、径流因子Q和风力因子W)的映射关系，即:

S=f(H,Q,W)

(1)

由于盐度与主要因子之间的关系十分复杂，咸度因子与其影响因子之间很难得出显示函数表达式. 目前，人工神经网络(artificial neural network，ANN)能够很好地处理这一类预测问题，并已应用于咸潮预测研究[22]. 本文拟采用反向传递人工神经网络(back propagation artificial neural networks，BP-ANN)，模拟取水口咸度与河口潮汐因子、上游径流因子和当地风力因子的关系，并通过上游来水过程预测取水口咸度.

在感潮河流供水安全保障规划中，当地水资源管理更加关注上游来水条件下取水时机问题. 根据国家饮用水源地水质标准要求，当含氯度超过250 mg/L被认为咸度超标，该阶段累积时间为超标时间. 因此，基于日尺度的咸度及其主要影响因子观测样本，通过BP-ANN模拟训练与验证，确定咸度预测模型结构和参数后，假定潮汐和风力因子维持原有定常条件，根据时尺度咸度观测样本和供水水质咸度标准限值要求，建立日均含氯度与日超标时间的映射关系，实现日超标时间与上游来水日均流量的连结，即:

Th=f(S)=f[f(Q|H,W)]

(2)

式中，Th为咸度日超标时间，S为日均含氯度，Q为日均流量.

通常情况下，日均含氯度值越高，每日咸度超标时间越长，属于正相关关系. 因此研究过程中，可采用了多种能够呈现正相关的曲线函数，拟合含氯度与超标时间之间的关系.

1.3 河库联合供水调度和供需平衡分析

受到咸度超标影响，感潮河流取水时机并不稳定. 取水口在咸潮影响时从河流避咸取水，每日避咸取水量首先供给水厂，若能满足区域需水，则多余取水量结合水库蓄水能力存储，若取水无法满足需水要求，则通过水库补充供给. 由于当地水库集水面积较小，枯水期供需平衡过程中不考虑水库集水区域的来水和蒸发损失等.

若河流取水口的日取水能力为Vc(m3),当地日需水量为Vd(m3),水库前一日蓄水量为W0(m3),水库总蓄水能力为Wc(m3),则取水口每日取水量Vw(m3)为:

Vw=Vc(1-Th/24)

(3)

若取水量不小于需水量，即Vw≥Vd,则实际日供水量Vs(m3)、水库日蓄水量VI(m3)和总蓄水量WT(m3)为:

Vs=Vd

(4)

VI=Vw-Vs

(5)

WT=W0+VI

(6)

式中，WT≤Wc. 若WT>Wc,则令WT=Wc.

若取水量小于需水量，即Vw

VO=Vd-Vw

(7)

WT=W0-VO

(8)

Vs=Vw+VO

(9)

式中，WT≥0. 若WT<0,则令WT=0,且有:

VO=W0

(10)

Vs=Vw+VO

(11)

式中，Vs

Ws=Vd-Vs

(12)

1.4 供水安全分析和临界流量

基于咸度预测和河库联合供水调度模块，可获悉某一频率或典型年上游来水条件下，避咸蓄淡过程的供需平衡状况. 从水资源规划和管理角度来说，探究满足未来发展需水要求的上游来水量即临界流量，成为了依赖感潮河段取水的城市供水安全关键问题之一. 一般来说，在水资源调配过程通常是基于供水保障率(包括供水水量保障率和供水天数保障率)分析供水安全[23]. 本文根据长系列上游来水过程，拟通过设定风险阈值和临界阈值，通过避咸蓄淡供水模型调度及供需平衡结果，分析上游来水的风险流量和临界流量.

有某年枯水期上游来水{Q0|Q0(i),i=1,2,…,N},N为总天数. 总需水为WD,通过咸度预测和河库联合供水调度之后，水库剩余蓄水量为WR(不含死库容)，总缺水量为Wa.设定风险阈值TR和临界阈值TC,有:

Wa≥0

(13)

WR/WD≤TR

(14)

Wa/WD≤TC

(15)

式(13)～(15)定量分析目的是，若存在缺水，则表明现有的河库联合供水系统对于该典型年上游来水而言存在供水风险；若不缺水，但联合供水调度期末的水库蓄水量较小时，考虑到咸度预测和来水过程等不确定性的影响，认为该典型年上游来水存在潜在的供水风险；若蓄水量和缺水量相对于需水而言在一个较小的范围之内，从供需平衡角度而言是可接受的供水风险，则该典型年上游来水被认为是期望上游来水的临界流量.

2 案例分析

2.1 区域概况和数据处理

如图2所示，粤港澳大湾区珠海及澳门位于珠江三角洲南部，其供水水源主要依赖于磨刀门水道取水. 作为珠江三角洲八大口门之一的入海通道，磨刀门水道属于典型感潮河流. 珠海市在磨刀门水道设有3个取水口，分别为广昌泵站、平岗泵站和竹洲头泵站，设计取水能力分别为100万、124万和80万m3/d. 3个泵站的取水不仅要保障珠海东城片区供水，而且还承担着对澳门供水的任务. 珠海市主要调咸水库为竹银水库，兴利库容为3811万m3；另在珠海市北部有若干个小型水库或塘坝构成了北部水库群，参与珠海东城片区和澳门的供水，合计兴利库容为2945万m3；总兴利库容为6756万m3.

图2 研究区地理位置及气象水文站、取水泵站、水库和供水对象

由于受到上游来水减少、河床下切和咸潮上溯等综合因素影响，磨刀门水道枯水期的咸度超标问题日益突出[24]. 尽管西江上游骨干水库群的压咸调度，能有效缓解河口地区咸潮上溯影响. 但是根据磨刀门水道主要取水口2010-2015年枯水期逐时咸度变化，广昌、平岗和竹洲头泵站取水口含氯度超标占比分别为70%、20%和11%，短时间尺度的咸度超标问题依然显著，如图3. 为此，珠海市加强了区域水库联合调蓄功能，基于磨刀门河流取水时机避咸供水，并充分挖掘当地水库调蓄能力取水蓄淡，形成了河库避咸蓄淡联合供水模式. 2018年，珠海市东部及澳门供水量约为3.61亿m3，枯水期需水量为1.81亿m3；参考当地区域发展及水资源规划，预测2035年需水量约为8.45亿m3，枯水期需水量为4.23亿m3.

图3 磨刀门水道主要取水口2010-2015年枯水期逐时咸度变化

磨刀门水道咸度观测开展较晚，考虑数据分析匹配，主要选用了2003年1月1日-3月31日、2003年10月1日-2004年3月31日、2004年12月1日-2005年3月31日和2005年10月1日-12月31日广昌泵站和平岗泵站的含氯度日均值，以及广昌泵站、平岗泵站和竹洲头泵站2011年10月1日-2012年3月31日的含氯度日均值和日超标时间. 上游径流因子以西江流域控制站高要和北江流域控制站石角的合计流量为输入，获取到1957-2008年(水文年)逐日流量过程. 潮汐因子以磨刀门水道入海口灯笼山潮位站的潮位过程为参照，风因子以最邻近研究区域中山气象站逐日风速风向为参照.

考虑到竹洲头泵站咸度数据较少，以广昌泵站和平岗泵站的逐日咸度因子含氯度日均值数据为样本，进行咸度预测模型训练和验证. 选取灯笼山日最低潮位、日平均潮位、日最大潮差、中山站风速风向、高要与石角的合计日流量，作为备选预测因子. 根据风对咸潮入侵影响的物理机制，以磨刀门水道流出河口方向为正建立直角坐标系，沿坐标系对风进行径向和切向分解，认为切向方向分量对咸潮有影响[21]. 潮位因子、流量因子和风速因子变化过程见图4.

图4 磨刀门水道咸度预测潮汐因子、径流因子和风因子变化过程

广昌泵站和平岗泵站备选预测因子与咸度之间相关性与互相关性见表1. 对比互相关程度，结合物理机制成因分析，流量、日平均潮位、日最大潮差、日最低潮位和切向风速对咸度影响滞后时间分别为+2、+1、+3、-2和+2 d. 除日平均潮位外，各预测因子与咸度的相关性均通过检验，且为负相关. 因此，确定磨刀门水道咸度预测模型的主要输入因子为流量、日最大潮差、日最低潮位和切向风速.

表1 磨刀门水道咸度和预测因子的相关性与互相关性

为了更好地了解当地需水增长过程中供需平衡以及上游来水过程差异对供水安全的影响，拟设定多需水情景下，分析不同来水系列的避咸蓄淡供水平衡状况. 基于现状用水(1.81亿m3)和未来规划需水(4.23亿m3)之间的差异水量，以其5%的水量等距增长，连续设定需水情景. 基于高要站和石角站1956-2008年枯水期(10月1日-3月31日)逐日合计流量，共获取52个来水系列过程. 根据枯水期平均流量系列进行统计模拟，获得每个来水系列过程流量均值特征及相应的设计频率，以此作为上游来水的特征流量. 在日供需平衡分析过程中，缺水情况下枯水期可供水量为需水量与缺水量之差的累积值，满足供需平衡情况下可供水量为枯水期累积供水量与末期水库结余量之和. 在风险流量和临界流量分析过程中，风险阈值和缺水阈值均设定5%. 若枯水期避咸蓄淡调度过程存在缺水或枯水期结束时水库蓄水量与总需水之比小于5%，则认为该来水过程流量均值为风险流量，其中枯水期总缺水量与总需水量之比小于5%及枯水期结束时水库蓄水量与总需水之比小于5%的来水过程流量均值为临界流量.

2.2 咸度模拟与预测

构建3层BP-ANN，将460组实测数据按照时间顺序分为模型训练期和验证期，分别有数据样本400和60组，模型拟合建立过程如下：

1)数据预处理. 考虑各因子数据类型差异，用Premnmx函数对输入输出数据归一化至[-1,1]后进行模拟，输出数据再进行反归一化处理.

2)模型构建及训练. 非线性传输函数为tansig，输出层采用线性函数purelin，选用traingdm算法，训练精度RMSE的阈值设定为[0.008,0.05][25]，优选出广昌泵站和平岗泵站模型结构分别为4-10-1和4-12-1. 训练后的模拟值与实际值相关程度较高，广昌泵站和平岗泵站的决定系数R2分别为0.95和0.84，构建的BP-ANN 模型收敛有效，如图5.

图5 磨刀门水道咸度预测的BP-ANN训练与验证

3)模型验证. 基于训练后模型结构和参数集，广昌泵站和平岗泵站在验证阶段的RMSE分别为0.46和0.11，决定系数R2分别能达到0.72和0.71，认为构建的BP-ANN是可接受的[26].

2.3 含氯度与超标实数关系拟合

由于超标时间是有上下边界的，其取值期间为0～24 h，拟合时基于边界约束对函数分段. 另外，磨刀门水道含氯度实测值变化范围跨多个数量级，拟合时对含氯度进行了对数变换. 广昌泵站和平岗泵站拟合结果如图6所示，中间曲线部分拟合的决定系数R2分别达到了0.88和0.90. 因此，可以用拟合曲线进行日均含氯度exp(x)与日超标时间y的转换. 拟合函数如下：

图6 广昌泵站(a)和平岗泵站(b)日均含氯度与超标时间的拟合曲线

广昌泵站为：

(16)

平岗泵站为：

(17)

由于竹洲头泵站缺少足够的咸度样本数据，无法采用BP-ANN进行模拟. 竹洲头站泵站与平岗泵站距离较近，为此建立竹洲头泵站与平岗泵站日均含氯度关系，如图7a所示，函数拟合的决定系数R2达到了0.96. 平岗泵站含氯度x转换为竹洲头泵站含氯度y的拟合函数为：

y=0.0001383x2+0.3267x+0.08153

(18)

基于竹洲头泵站咸度观测数据，拟合含氯度exp(x)与超标时间y的关系，结果如图7b，中间曲线部分拟合的决定系数R2达到了0.92. 拟合函数为：

图7 平岗泵站与竹洲头泵站日均含氯度关系(a)及竹洲头泵站日均含氯度与超标时间(b)的拟合曲线

(19)

2.4 不同需水情景下及典型年的供需平衡

基于避咸蓄淡供水模型，不同需水情景下所有年份枯水期来水过程的可供水量如图8a. 在需水量为1.81亿～2.54亿m3时，所有来水系列可供水量基本一致，不需要水库参与调蓄供水，3个取水口避咸取水直接供水就能满足此阶段需水要求. 之后的需水增长，对于部分年份来水过程来说需要水库参与调蓄供水. 当需水增长到情景3.75亿m3时，1964和1998年来水过程的可供水量小于需水量，开始出现了缺水. 之后随着需水继续增长，出现缺水的来水年份数也不断增加. 到2035年需水4.23亿m3时，共有21 a的来水过程会出现缺水. 因此，需水量3.22亿m3为当地河流取水直接供水满足需水要求的临界需水，可视为枯水期资源性缺水的临界需水量；需水量为3.75亿m3为当地避咸蓄淡联合供水满足需水要求的临界需水，可视为枯水期工程性缺水的临界需水量. 在水库参与调蓄供水过程中，随着需水量增加，可供水量也逐渐增加. 但可供水量增长速率是有差异，且部分年份的增长速率相继变缓，即表明对于该来水过程而言存在最大可供水量的临界需水. 实际上，由于受到当地河库联合调蓄能力的限制，如果未来需水持续增加，对于任何来水过程都存在最大可供水量的临界需水.

不同需水情景下长系列来水过程的供需平衡结果如图8b所示. 在需水量小于3.75亿m3时，所有年份来水过程均能满足需水要求. 当需水量由3.75亿m3增加到4.23亿m3时，出现缺水的时间由2 a增加到21 a. 到2035年枯水期需水达到4.23亿m3，即使采取避咸蓄淡联合供水模式，依然有近40%的来水过程无法保障当地的需水要求.

图8 不同需水情景下供需平衡：供水与需水(a)；供水与缺水(b)

选择上游来水枯水期平均流量频率接近为P=90%的典型年，包括1962年(P=88.7%)和2004年(P=90.6%)，分析2035年枯水期需水下供需平衡过程，超标时间和供需平衡情况见图9和图10. 1962年，10-11月及3月中下旬的流量较大，12月中旬至3月上旬持续维持较小流量. 2004年，除了3月下旬流量较大外，其他时间段流量比较均匀. 1962和2004年，枯水期缺水量分别为2681万和1297万m3，缺水时间分别为42和27 d. 尽管2个典型年的枯水期平均流量差别较小，但是缺水量和缺水时间差别显著，表明来水过程特征对避咸蓄淡供水结果影响较大. 由于预设枯水期之前水库处于蓄满状态，即使枯水期前期上游来水流量较小，咸潮影响下日取水时数较短，通过水库供水能满足一定时段内的供水要求. 但是较小流量在中后期维持时间较长时，水库水量耗竭且取水供水后的剩余量补库不足，则容易出现缺水现象.

图9 典型年枯水期流量及超标时间：1962年(a)；2004年(b)

图10 典型年枯水期供需平衡结果：1962年(a)；2004年(b)

2.5 风险流量和临界流量分析

不同需水情景下风险流量及临界流量见图11a. 由于从安全的角度考虑了枯水期末水库最小剩余蓄水量的要求，需水量在3.50亿～4.23亿m3时，供水存在风险的时间分别有2、4、5、13、18、21和26 a，相比于缺水时间分别增加了2、4、3、9、11、5和5 a. 需水量为3.50亿m3和3.63亿m3时，尽管不存在缺水年份，但有2 a由于枯水期末水库剩余蓄水量较小，认为存在供水风险. 基于风险阈值和缺水阈值识别出临界流量年数，分别占供水风险年数的3.8%、7.6%、9.6%、25.0%、34.6%、40.4%和50%. 在2035年需水量为4.23亿m3时，枯水期平均流量最大的年份也存在供水风险，进一步说明不是枯水期平均流量越大供水越安全，还与枯水期流量分布有关.

不同需水情景下临界流量统计特征如图11b和表2所示. 整体上来说，随着需水增长，对上游来水期望的临界流量越大. 大约分为3个阶段：枯水期需水量为3.50亿～3.62亿m3时，临界流量均值在2373～2411 m3/s之间，对应设计频率为96.8%～98.1%，相对变幅为-12.9%～12.9%；当需水量增至3.74亿～3.98亿m3时，临界流量均值在2723～2828 m3/s之间，对应的设计频率范围为76.7%～82.3%，相对变幅为-11.6%～15.7%；当需水量增至4.10亿～4.23亿m3)时，临界流量均值为3322～3372 m3/s，对应设计频率为50.0%～52.2%，相对变幅为-13.4%～13.6%.

图11 不同需水情景下风险流量(a)与临界流量(b)

表2 不同需水情景下临界流量统计特征

3 结论

1)BP-ANN能够较好地模拟感潮河流的咸度与上游径流因子、河口潮因子和当地风因子的映射关系，并结合取水口含氯度与超标时间的曲线函数模拟，能够建立上游来水-咸度-超标时间之间的联系. 结合当地水库调蓄能力和需水要求，可得到枯水期上游来水过程下的当地取水口取水时数、取水量以及水库调蓄和区域供需平衡状况. 通过设定供水风险阈值和临界阈值，基于长系列枯水期来水过程识别出了不同需水规模的风险流量和临界流量. 因此，通过构建耦合了咸度预测、河库联合供水调度和供水安全分析模块的避咸蓄淡模型，为依赖感潮河流为水源地的区域供水安全分析提供了一种整体思路和决策方法.

2)对于依赖珠江三角洲磨刀门水道取水的粤港澳大湾区珠海东部及澳门供水来说，资源性需水临界点和工程性需水临界点分别为3.22亿和3.75亿m3. 当2035年规划需水量4.23亿m3时，上游枯水期平均流量频率约为90%的典型年(1962年)实际来水过程下的缺水量为2681万m3. 当枯水期需水量由3.50亿m3增至2035年4.23亿m3时，上游临界流量均值由2373 m3/s上升到3372 m3/s，对应设计频率由98.1%变化为50.0%.

3)整体上来说，需水规模越大，对上游来水期望的临界流量越大. 西江上游枯水期调度对于保障珠江河口地区供水安全的作用是显著的，但短时间尺度咸潮超标问题依然存在，需要充分挖掘河口地区水库供水潜力，通过避咸蓄淡河库联合供水才能更好地保证2035年需水要求. 某一需水情景下的临界流量也不是唯一的，还与枯水期流量分布有关. 如何分解来水过程的流量分布特征以及河流取水能力及水库调蓄能力变化，对感潮河流区域供水安全的影响，有待于进一步深入研究.