1. 引言
巢湖是我国五大淡水湖和“三河三湖”之一，具有防洪、灌溉、供水、航运、水产、旅游等多种功能，同时水多、水少、水滞、水脏等问题交织并存，新老水问题极其复杂 [1] [2] [3] 。巢湖流域形状呈东西长、南北短，发源于大别山余脉和江淮分水岭，最大河流杭埠河自西向东流入巢湖，另有南淝河、派河、十五里河、白石天河、柘皋河和兆河等河流呈放射状注入，并由裕溪河、牛屯河排入长江。巢湖闸以下裕溪河水系主要支流有清溪河、西河等，西河、兆河、裕溪河也是巢湖重要的引江水道。巢湖流域现有水文(位)站20余个，湖区入汇的杭埠河、南淝河等较大支流上设立水文站较早，具有较长的实测流量资料；派河、十五里河、白石天河、兆河等较小河道流量资料较为缺乏(未设站，或设站时间较短，或仅观测水位)。
巢湖入湖河流的径流特征及变化分析是开展水旱灾害治理、水环境治理、蓝藻水华防控、水生态修复的重要基础，也是分析河湖污染通量、研究水质水量响应关系、开展“引江济淮”工程的优化调度、提升巢湖综合治理水平的重要依据。长期以来，因主要聚焦防洪减灾及汛期洪量，巢湖缺乏入湖径流系列数据，目前入湖39条河道仅有2个系列不完整的水量监测断面，制约水多、水少、水滞、水脏等巢湖综合治理规划设计和河长制管理工作。为掌握入湖径流变化规律，支撑巢湖综合治理，提升巢湖科学调度水平，开展巢湖主要入湖河道长系列逐月逐旬径流回溯模拟研究不仅是必要的，也是急迫的，具有重大的应用价值和现实意义。
2. 基于水量平衡方法反演巢湖总入湖径流过程
2.1. 反演方法
巢湖流域的入湖河道和巢湖湖区均可视为一个系统，巢湖入湖河道径流汇入湖区，湖区同时承受降雨增益和蒸发损失，经过湖区调蓄，汇入巢湖闸的径流为输出结果。因此，巢湖流域的水量平衡计算可以利用多输入单输出(Multi-input and single-output, MISO)系统模型进行分析。基于MISO模型的基本原理，建立巢湖流域水量平衡的MISO模型如图1所示，主要选取了巢湖控制流域面积较大的主要河道，将其余小河道概化为河道下游圩区和巢湖周边区间，并进一步考虑巢湖周边取用水。
巢湖闸设站后，具有长系列的出闸流量过程；兆河闸设在巢湖的入湖河道兆河上，也具有长系列的流量过程，由兆河汇入巢湖的流量被记为正值，由巢湖汇入兆河被记为负值。采用巢湖闸和兆河闸整理后的1967~2021年共55年的日流量资料作为巢湖的总出湖流量资料；采用巢湖的4个水位站资料均化后的日水位数据(采用巢湖闸上游水位、忠庙水位、塘西水位、槐林水位)，依据巢湖的高程(H)-面积(A)-容积(V)曲线，从而反演巢湖的库容变化。本文搜集整理了巢湖流域7个国家基本气象台站，涵盖1955~2021年共67年的气象资料，主要有日降水量、日平均气温、日最高气温、日最低气温和日蒸发量，其中蒸发量分别有大型蒸发皿和小型蒸发皿的观测值。本文还收集到174个自动气象观测站，涵盖2015~2021年共7年的气象资料，包括日降水量、日平均气温、日最高气温和日最低气温，无蒸发观测资料。国家基本气象台站中仅有巢湖站位于巢湖湖区，故采用该站的长系列日蒸发资料作为湖区的面蒸发系列，并基于泰森多边形方法推求巢湖湖区的面降水系列。
Figure 1. Diagram of the MISO model in evaluating water balance in Chao Lake basin
图1. 巢湖流域水量平衡演算的MSIO模型概化图
2.2. 巢湖总入湖径流过程及变化特征分析
首先研究不同历时状态下的巢湖国家基本气象台站观测降水量与湖区面降水量回归关系，发现国家基本气象台站的降水观测值与自动站推求的湖区面降水系列呈现显著的正相关性，且显著性水平p值均低于0.05。同时，可以发现，月降水量的相关特征最优，日降水的相关性最差，但是线性回归方程仍然取得了较好的拟合效果，R2值超过0.86。鉴于线性回归方程良好的拟合效果，将该回归关系应用到1955~2014年，基于巢湖国家基本气象台站的实测降水系列，反演得到巢湖湖区的长系列面降水过程，记为方案一。仅采用国家基本气象台站表征湖区降水系列，记为方案二。基于该湖区面降水、蒸发、水位及湖区周边取用水资料，采用MSIO模型反演得到了两种方案下巢湖的总入湖径流系列。
采用Mann-Kendall (M-K)趋势检验、Petti突变检验和周期性分析方法分析巢湖年入湖径流的变化特征，发现巢湖总入湖年径流过程呈现先上升后下降的特点，但是在0.05显著性水平下变化不显著；巢湖的年径流量系列存在突变，而且其突变时间在1975~1980年左右；在研究时段，巢湖湖区的年降水量在研究时段内没有显著性周期变化，但在其中某些时段存在一定的周期，例如在2000年左右存在显著的周期性，如图2所示。
Figure 2. Periodic analysis of annual runoff series in Chao Lake
图2. 巢湖总入湖年径流系列的周期性分析结果
3. 基于卫星遥感反演巢湖流域长系列气象数据集
巢湖流域的国家基本气象台站，涵盖1955~2021年资料，但是空间密度低；流域内近年来设立了174个自动气象观测站，但是仅观测2015~2021年的日气象系列。因此，本文尝试融合有限的地面气象观测数据、国际上较新的长系列高精度Multi-Source Weighted-Ensemble Precipitation Version 2 (MSWEP-V2)卫星集成降水数据集和欧洲中期天气预报中心的ERA5气温数据，首先通过基于分位数映射的日偏差校正(DBC)、基于月尺度的回归校正(LRBC)和等率校正(RBC)等3种偏差校正方法，对遥测栅格降水和再分析气温日系列进行校正，再采用季节性贝叶斯模式平均(Bayesian Model Averaging，BMA)方法求解各偏差校正系列的后验分布优选相应权重，得到融合多种偏差校正模式的日降水、气温系列，最后选取2个试验子流域分别基于新安江、GR4J和HMETS模型评估其水文模拟适用性及精度，关于MSWEP-V2和ERA5数据集的详细信息及具体校正方法，详见文献 [4] 。主要结论如下：
1) MSWEP V2降水和ERA5气温数据集在巢湖流域均存在一定的模拟偏差，故在流域水文模拟应用领域有必要开展气象数据偏差校正；偏差校正模型需要可靠的观测资料进行率定，该流域具有空间密度较高的自动气象站，为构建校正模型提供了良好的适用条件；对于地面观测站点稀少的地区，降水受地形、气候等多种因素影响较大，其偏差校正方法的稳定性有待进一步探讨。
2) 4种偏差校正方法对气温变量进行校正后，整体上均得到较低的模拟误差，其中DBC方法对气温极值的模拟能力最优，BMA方法次之；对于均值气温，DBC、RBC和BMA方法均能得到优良的校正效果，校正后各月均温与实测系列的差异一般在0.5℃以内，日气温的相关性系数均超过0.9。
3) 日降水相较于气温系统误差偏大，LRBC方法对降水的校正能力相对较差，而BMA、DBC和RBC方法均取得较好的校正效果，这3种校正模式下大多数评价指标的相对偏差低于5%；基于国家基本气象台站长系列资料的泰勒图和回归分析结果显示，BMA方法优于其他模式，日降水的相关性系数接近0.8，月降水量的相关性系数超过0.94。
4) 采用BMA模型偏差校正后的气象数据系列，分别驱动新安江、GR4J和HMETS水文模型，结果显示在2个试验子流域均取得较好的模拟效果，率定期和检验期的Kling-Gupta efficiency (KGE)系数超过0.67，纳西效率系数(NSE)和相对偏差(RB)等指标显示水文模型的水量误差较小。总体而言，BMA方法校正的日降水、气温系列能够满足水文模拟的应用要求。
4. 桃溪水文站日径流系列模拟与插补
巢湖流域仅有两个系列不完整的水量监测断面，其中较为关键的丰乐河桃溪水文站也存在较多的缺测漏测现象，而巢湖流域反演的长系列降水和气温数据集具有时空分辨率高、数据质量好的特点，故可以通过丰乐河气象水文资料建立流域水文模型，实现水文资料的插补延长和资料稀缺流域的径流重塑模拟。杭埠河是巢湖流域最大的入湖河道，在杭埠河干流中上游建有龙河口水库，丰乐河是杭埠河最大的支流，而且其与杭埠河干流的交汇点接近巢湖湖区，桃溪站是丰乐河中上游的水文站。因此，通过重构桃溪水文站的径流系列不仅可以回溯模拟丰乐河的径流过程，还能为其他缺资料流域水文模拟过程的参数移植提供参考。本文采用新安江模型和GR4J模型率定丰乐河流域的径流过程，并采用水文模型插补桃溪水文站的缺测系列。对于水文模型的参数率定方法，具体采用复合型混合演化(SCE-UA)算法优选水文模型参数，该算法是一种全局优化算法，它集成了随机搜索算法、单纯形法、聚类分析及生物竞争演化等方法的优点，能有效处理目标函数反映面的不敏感和不凸起等问题，且不受局部最优点的干扰 [5] [6] 。
桃溪水文站控制流域的水文模型精度普遍较好，其中GR4J模型率定期的KGE系数一般超过0.85，模拟的旬径流和月径流与实测系列的相关性系数一般超过0.72；新安江模型也取得了良好的模拟效果，率定期的KGE系数一般超过0.7。图3给出了桃溪水文站控制流域GR4J模型日尺度、旬尺度和月尺度的模拟与实测径流回归分析图，从图中可以看到模型的日径流模拟系列与实测系列较为接近，在旬尺度和月尺度均得到较高的相关性系数，仅在第III子期较低(但是旬和月尺度均超过0.69)，而在其他三个子期日、旬和月尺度的相关性系数一般在0.7~0.8之间。
总体来看，丰乐河流域建立的流域水文模拟总体上精度较高，能够满足巢湖流域无资料地区入湖河道的径流回溯模拟要求，故采用桃溪水文站控制流域的降水和气温数据驱动新安江模型和GR4J模型，分别模拟得到两套日径流系列；考虑到两个水文模型的模拟效果均较好，故最后采用两套日径流系列的均值替补缺测系列，得到桃溪水文站径流过程。
5. 基于产汇流经验计算模拟巢湖天然入湖径流系列
把巢湖流域分为山丘区和圩区，各分块山丘区通过次降雨径流关系得到净雨过程，再通过净雨与模型单位线复合得出各分块的径流过程线，再与圩区排水过程逐日叠加，得巢湖闸上逐日径流过程线。山丘区径流的计算过程分为产流和汇流两部分，产流和汇流计算过程如下。全流域目前雨量站分布较多，但是部分雨量站建设较晚，不具有长系列的降雨资料，部分雨量站监测数据不连续，没有完整的长系列资料。本报告采用两种方案选用降水数据：方案I为第3节推求的1979~2021年栅格化降水数据集，通过选取各子流域范围内所有栅格的降水进行均化，得到各河道的面雨量过程；采用的方案II为雨量站，根据雨量站的位置和资料的完整性进行选择，每条河道上、中、下游均会选择合适的雨量站，所选的雨量站尽量分布均匀，本次共选择雨量站28处，根据泰森多边形计算各雨量站的权重。利用所选用的雨量站计算各分块面积上的场次降雨，采用次降雨径流关系计算场次径流，再分配到每一天，得到逐日的流量过程，次降雨径流关系采用桃溪站和合肥站的次降雨径流关系曲线。以每条河道作为一个分块，巢湖闸上分为8块。山丘区模型单位线，采用无资料地区的模型线，这种模型是在特定前流域特性(K)和时段3b (Δt = 3 h)，降雨强度为50 mm，流域面积为1000 km2的情况下做出的。经过换算做出了不同分块面积上，不同流域特性(K)，时段为1 d (Δt = 24 h)的模型线。
Figure 3. Regression results of GR4J-simulated and observed runoff series in control basin of Taoxi hydrological station
图3. GR4J模型在桃溪水文站控制流域模拟与实测径流的回归分析结果
Figure 4. Natural total inflow derived by yielding and routing method and observed annual runoff in Chao Lake
图4. 基于产汇流计算得到的巢湖天然总入湖与实际入湖年径流
基于上述方法，采用栅格降水数据集和雨量站数据分别模拟得到巢湖各入湖河道的天然径流系列后，通过与圩区排水系列叠加，得到天然状态下的巢湖总入湖系列(见图4)。从该图看出，采用雨量站点数据和栅格降水得到的总入湖径流系列非常接近；由于推求栅格化降水数据时采用了巢湖流域所有搜集到的实测雨量数据进行校正，各站点月尺度的栅格降水与实测降水的相关性系数达到0.8，校正过程还使用年降水总量进行了质量控制，故经过流域面平均后得到的年系列与实测数据差异进一步缩小，这一结果也同时验证了第3节反演栅格降水数据集的合理性。为了对比分析天然径流系列的合理性，图4进一步给出了基于产汇流计算的天然总入湖与水量平衡反演实际总入湖系列的逐年径流量，可以看出天然入湖与实际入湖系列的年际变化趋势整体上接近，且能较好地刻画丰水年、枯水年和平水年。从图中进一步看出，天然径流系列与实测径流系列存在一定差异，这是由于巢湖流域水库调蓄、工农业取用水等影响，故有必要进一步考虑人类活动影响、推求巢湖各河道的入湖径流系列。
6. 巢湖入湖河道径流方案比选与基于人类活动影响的校正优化
6.1. 多方案比选分析
采用常规水文比拟法和考虑降水修正的水文比拟法分别推求南淝河、白石天河、派河、柘皋河和十五里河等五条集水面积超过100 km2的入湖河道长系列径流过程，并对两种方案下的长系列入湖径流过程进行集成；圩区排水过程大圩按照0.6~1.1 m3/s/km2进行计算；最后通过水量平衡进行校正，从而得到优化后的各河道入湖径流系列。由于巢湖流域39条入湖河道中仅有7条河流集水面积超过100 km2，故将剩余1259 km2概化为其他区间。
经过比较分析，常规水文比拟法得到的径流系列在南淝河和白石天河等河道难以满足降水–径流关系，这是由于各河道的下垫面条件差异较大，难以直接采用面积缩放方式获得径流过程，故考虑降水因素对模拟入湖河道径流系列具有重要意义。进一步对比考虑降水修正的水文比拟法与产汇流计算的天然径流过程，发现杭埠河的年际变化趋势基本一致，其他河道的差异较大，如图5所示。从该图看出，尽管两种方案得到的杭埠河年径流年际变化趋势接近，但是各年份的年径流总量存在差异。
Figure 5. Comparison of annual runoff during 1968~2021 based on two different methods in Hangbu River
图5. 两种方案下杭埠河1968~2021年径流的对比结果
通过综合分析和多方案比选，发现各方案各有优劣，水文比拟法和降水修正法的计算量较小，对有资料地区的河道径流模拟精度较高，例如能够充分利用杭埠河的观测系列，但是由于巢湖各入湖河道的产汇流条件存在显著差异，这两种方案在其他河道的应用效果较差，该方案推算的总入湖与实测系列差异较大，难以直接应用于回溯模拟总入湖径流系列。采用产汇流计算得到的天然入湖系列能够较好地捕捉总入湖径流过程，各河道及圩区叠加的天然入湖量与实测总入湖的年际变化趋势接近，但是未考虑水利工程和工农业用水影响，对部分河道的模拟结果存在偏差。
尽管上文推求了天然状态下的各河道入湖径流，为了适应现状条件下支撑推求入湖河道污染通量的要求，以满足调度应用和水资源管理等实际需求，当前迫切需考虑人类活动的影响量。综合上述分析结果，本文建议以经验产汇流计算的天然入湖径流系列为基准，进一步考虑流域外引调水工程、湖区周边工农业取用水等人类活动，给出优化校正后的巢湖入湖河道径流系列。
6.2. 考虑外流域调水影响校正净雨量和径流系列
把巢湖流域分为山丘区和圩区，各分块山丘区通过次降雨径流关系得到天然净雨过程，再通过净雨与模型单位线复合得出各分块的径流过程线，再与圩区排水过程及水库下泄过程逐日叠加，得巢湖闸上天然逐日径流过程线。由于防洪排涝工程的建设，流域防洪排涝能力不断提高，根据各典型年实际统计相应的排涝能力，初步统计，现状圩区泵站排水模数一般在0.6~1.1 m3/s/km2。
基于两套降水数据集，通过以上方法进行产汇流计算，得到各分块山丘区的逐日径流过程，各分块径流过程叠加圩区排水过程、大型水库泄水过程和湖面的来水过程，得到巢湖闸上逐日径流过程。考虑到兆河的河道较为平缓，且来水交错复杂；以兆河闸为界，兆河闸下83 km2直接入巢湖，兆河闸上537 km2入西河，难以直接反演。本文以柘皋河为参照流域，首先采用面积比法缩放，折算后得到兆河闸上的径流系列，与兆河闸的汇入流量叠加，得到兆河的入湖径流系列。两种降水方案下得到的校正后径流系列，分别记为方案I和方案II。两种降水方案得到的入湖径流基本接近(均基于实测数据进行了校正处理)，主要在洪水年和枯水年的若干月份存在差异。例如，对于杭埠河而言，不同方案的径流系列在1991年和2020年存在一定不显著的差异，在其他年份普遍非常接近。整体来看，两种方案得到的径流结果较为接近，方案I采用栅格降水系列得到面降雨进而推求入湖径流能更充分地考虑气象事件的空间异质性，方案II采用雨量站资料与以往的设计规划报告较为接近。进一步比较两种方案的年际变化，发现方案I对捕捉特大/小年份的径流过程较为有优势，这是由于栅格降水数据集能较好反映巢湖流域的空间异质性及降水事件的时空过程。基于以上分析结果，对于1979~2021年系列，本文采用栅格降水数据(方案I)得到的径流过程；对于1951~1978年系列，由于暂无栅格降水数据，采用雨量站(方案II)推求的径流过程。图6给出了上述方法得到的杭埠河年径流的Pettitt突变检验结果，可以看到回溯模拟的杭埠河年径流系列在2010年存在显著突变。
Figure 6. Abrupt point detection analysis of Hangbu River annual inflow runoff
图6. 杭埠河入湖年径流的突变检验分析结果
6.3. 合理性分析
为进行成果合理性说明，首先比较本文推求的径流回溯模拟过程和水量平衡方法得到的实测过程，差别最大的年份为1969、1991和2020等大水年，这几个大水年7月份均发生了破圩，进一步分析发现月平均差值最大的月份也为7月份，故结果基本合理。根据《安徽省第三次水资源调查评价》成果，巢湖流域1956~2016年多年平均地表水资源量为59.79亿m3，50%、80%和95%保证率下的地表水资源量分别为66.39亿m3、35.73亿m3和23.65亿m3。杭埠河流域多年平均地表水资源量为21.91亿m3，50%、80%和95%保证率下的地表水资源量分别为24.84亿m3、12.17亿m3和8.34亿m3。本文推求的巢湖1951~2021年多年平均入湖径流量为43.18亿m3，量级相比水资源调查评价的结果偏小，这是由于本文仅考虑巢湖闸上游的入湖径流量，未考虑巢湖湖区、闸下游及湖区的径流量，故偏小结果是合理的。进一步对比杭埠河的多年平均径流量，发现本文推求的径流量为20.82亿m3，1956~2016年径流量为20.98亿m3，与水资源调查评价的结果略微偏小；但是《安徽省第三次水资源调查评价》成果给出的地表水资源量是指河流、湖泊、冰川等地表水体中由当地降水形成的、可以逐年更新的动态水量，是天然河川径流量，本文推求的杭埠河径流系列考虑了水库调蓄影响和工农业耗水过程，结果适当偏小是合理的。
7. 巢湖总入湖年径流频率分析
考虑到不同水文频率分布线型和参数估计方法可能对计算结果有较大影响，本项目同时采用P-III/CF和P-III/LM两种不同的频率分析模型开展巢湖总入湖年径流频率分析计算，相关计算成果见图7和图8。其中，考虑国内通常做法和地区综合需要，对P-III/CF采取Cs/Cv = 2。关于两种频率分析模型的介绍，参见文献 [7] [8] 。结合图7看出，采用中国规范推荐的P-III/CF模型时，巢湖总入湖年径流1%设计频率下的径流量为103.7亿m3，2%设计频率下的年径流量为94.2亿m3，10%设计频率下的年径流量为70.4亿m3，20%设计频率下的年径
Figure 7. Hydrological frequency results of annual inflow runoff of Chao Lake by P-III/CF model
图7. 巢湖总入湖年径流P-III/CF模型频率曲线结果
Figure 8. Hydrological frequency results of annual inflow runoff of Chao Lake by P-III/LM model
图8. 巢湖总入湖年径流P-III/LM模型频率曲线结果
流量为58.7亿m3；若采用线性矩法估计P-III分布函数，设计年径流量普遍减少6%~16%，例如1%设计频率下的年径流为88.7亿m3，减少了14.5%。
本文进一步开展了各入湖河道的年径流频率分析，发现采用中国规范推荐的P-III/CF模型时，杭埠河入湖系列年径流1%设计频率下的径流量为48.5亿m3，2%设计频率下的年径流量为44.2亿m3，10%设计频率下的年径流量为33.4亿m3；若采用线性矩法估计P-III分布函数，设计年径流量普遍减少3%~7%。
8. 结论
巢湖入湖河道的长系列径流回溯模拟结果对支撑“引江济淮”工程具有重要意义，本文首先采用水量平衡方法推求了巢湖总入湖径流系列，进一步融合卫星遥测数据、气候模式再分析数据和地面观测数据，分别采用DBC、LRBC、RBC和BMA等四种偏差校正方法对卫星遥测数据进行校正，得到了巢湖流域1979~2021年的高分辨率气象数据集，并验证了其各项精度满足要求。采用新安江模型和GR4J模型对桃溪水文站控制流域开展了模型率定，并验证了模拟精度满足要求，各模型的参数在合理区间，具有物理意义，利用该模型插补了桃溪水文站的日径流系列。基于卫星遥感反演及雨量站数据，从产汇流过程推求巢湖总入湖及各河道的天然入湖径流系列，把巢湖流域分为山丘区和圩区，各分块山丘区通过次降雨径流关系计算净雨过程，考虑人类活动影响校正净雨量，再通过净雨与模型单位线复合得出各分块的径流过程线，再与圩区排水过程、大型水库泄水过程和湖面的来水过程逐日叠加，得巢湖闸上逐日径流过程线，并验证了该结果的合理性。
考虑引调水工程、圩区排水过程、工农业耗用水、污水厂排污等，基于多种方案推求了巢湖的总入湖及各入湖河道的径流系列，通过多方案比较和优化，给出了各河道的1951~2021年入湖径流过程，并区分了天然径流量和人类活动影响后的径流系列。经过回溯模拟，巢湖的年入湖径流量约为43.18亿m3；采用P-III/CF和P-III/LM两种不同的频率分析模型，开展了年径流频率分析计算工作，采用中国规范推荐的P-III/CF模型时，巢湖总入湖年径流1%设计频率下的径流量为103.7亿m3，2%设计频率下的年径流量为94.2亿m3。
基金项目
巢湖生物资源调查及生态修复示范工程(基础调查研究部分) (编号：2021BFFBZ01571)、国家自然科学基金(编号：52009091；52242904)、中央高校基本科研业务费专项项目(编号：2042022kf1221)、武汉大学研究生导师育人方式创新项目、武汉大学学位与研究生教育教学改革研究项目、武汉大学“教育教学改革”建设引导专项(本科教育质量建设综合改革子项目)、武汉大学研究生研究学分课程建设项目等资助。