摘要
为解决直接利用原始DEM提取的模拟河网与实际水系存在偏差、影响子流域划分的问题,论文基于区域实际影像中的河流位置,提取研究流域的主河道;依据主河道栅格高程与垂直主河道方向的洼地栅格高程之差的最大值,降低主河道所在DEM栅格的高程值,并采用河网密度法确定河源的最佳集水面积阈值。结合我国一级小流域面积范围标准(50~300 km2),以及每个子流域出口位于河网干流上的准则,确定合适的子流域划分数量及面积,并对不同的子流域划分方案下分布式新安江模型的日模和次模结果进行比较。桃溪流域应用结果表明:将实际水系主河道所在DEM栅格的高程降低40 m,基于河网密度法确定河源的最佳集水面积阈值为0.81 km2后,所提取的模拟河网更加符合实际。子流域最佳集水面积阈值为108 km2,划分的子流域数量为7个,子流域的平均面积为213.2 km2。对主河道所在DEM栅格的高程降低处理后划分的7和11个子流域分别记为方案A和B,未对主河道所在DEM栅格的高程降低处理划分的7个子流域记为方案C。利用桃溪流域1982—2022年水文气象资料对日径流过程进行模拟(1982—2012年为训练样本,其余为验证样本),方案A、B、C在训练期的确定性系数分别为0.87、0.87、0.86,验证期分别为0.90、0.91、0.90,模拟精度基本一致。选取20场洪水进行次洪模拟,以洪峰、洪量和峰现时间误差为精度评定指标,方案A、B、C下的合格场次分别为19、19、16场;以确定性系数为精度评定指标,方案A、B、C下的平均确定性系数分别为0.93、0.94、0.93。方案A和B的洪水模拟合格率高于方案C,而划分为7个子流域的方案A相较于11个子流域的方案B在精度上无明显差异。可见对主河道所在DEM栅格进行高程降低处理,以子流域面积处于我国一级小流域面积范围标准(50~300 km2)内、子流域出口位于干流上的准则进行子流域划分是合理可行的。
Abstract
To address deviations between the DEM-derived simulated river networks and actual hydrological systems that impacted sub-basin delineation, this study first extracted primary river channels from the study area using geospatial imagery. Then we adjusted elevation values in the original DEM by vertically reducing river grid elevations along the identified channels. The reduction magnitude was determined by the maximum elevation difference observed between the main river grids and the adjacent depression grids perpendicular to the fluvial pathways. The optimal watershed area threshold for the river source was determined by the river network density method. Subsequently, combined the standard range of China's first-level small watershed division (50-300 km2) and the criterion that each sub-watershed (SW) outlet was located on the mainstream of the river network, the appropriate number and SW area were determined. Finally, by comparing the daily and hourly results of the Xin'anjiang model distributed in different SW division schemes, the rationality of the SW division was verified. The application results in Taoxi watershed showed that: Lowering the elevation of the DEM grid where the actual river system's main rivers pass by 40 m and determining the optimal watershed area threshold for the river source by the river network density method, the simulated river network extracted was more consistent with the actual water system. The optimal area threshold for SWs was 108 km2, and the number of SW was 7, with an average area of 213.2 km2. The 7 and 11 sub-basins divided by the DEM grid elevation reduction treatment after the main river course was passed were scheme A and B, respectively, while the 7 sub-basins divided without the DEM grid elevation treatment after the main river course was passed were scheme C. The daily discharge process was simulated using the hydro-meteorological data of Taoxi watershed from 1982 to 2022 (the years from 1982 to 2012 were used as the training samples, and the rest were used as validation samples). The determination coefficients of the training samples of schemes A, B and C were 0.87, 0.87 and 0.86, and the validation coefficients were 0.90, 0.91 and 0.90, respectively. The simulation accuracy was basically consistent. A total of 20 floods were selected for secondary flood simulation, and the flood peak, flood volume and error of the time of flood peak were used as precision evaluation indexes. The qualified floods under scheme A, B and C were 19, 19 and 16, respectively. Using the certainty coefficient as the accuracy evaluation index, the average certainty coefficient of the next flood simulation in scheme A, B and C were 0.93, 0.94 and 0.93, respectively. Scheme A and scheme B had higher qualified rate of flood simulation than scheme C, and scheme A divided into 7 sub-basins had no obvious difference in accuracy compared with scheme B divided into 11 sub-basins. It can be seen that it is reasonable and feasible to reduce the elevation of the DEM grid through the main river channel and divide the SW according to the criteria that the SW outlet is located on the main river channel and the SW area is within the classification range of China's first-class small watershed (50-300 km2).
数字高程模型数据DEM[1]与GIS技术相结合,可识别更加贴近实际的流域水系特征。栅格型DEM由于处理起来简便有效,被广泛应用于流域水文过程模拟。对于分布式水文模型而言,精确获取流域特征信息,特别是河源位置的准确判定以及子流域的合理划分,对模型最终模拟结果的准确性和可靠性具有至关重要的影响。通过综合应用实际遥感影像数据与DEM数据,优化中小流域的子流域划分方法,可提升水文过程模拟精度,为水资源管理与洪水预警等提供科学依据。
众多学者在基于DEM提取河网方面开展了大量的研究:Yang等[2]提出了一种基于5 m分辨率DEM的半自动沟谷网络提取方法,该方法通过迭代沟谷汇流阈值深化沟谷网络,并结合“相对沟深(RGD)”的地形测量方法,可在黄土高原丘陵区提取较高精度的沟谷网络;卢庆辉等[3]提出了一种结合Priority-Flood算法漫水思想和改进D8算法的流向算法,有效地解决了平地区域生成平行河网的问题;巫晓玲等[4]基于5 m分辨率DEM,提出了一种基于沟头点群汇流分布的沟谷汇流阈值确定方法,结果表明在黄土高原不同地貌类型区适应性较好,但仍可能存在生成平行河网的问题。而根据实际影像确定河流位置,结合DEM提取河网和子流域的研究尚不多见。河网密度法在确定模拟河网提取阈值的研究中已有丰富成果:李照会等[5]分别在贡曲流域、辰清河流域和藤条江流域拟合了河网密度与集水面积阈值、河源密度与集水面积阈值之间的关系,结果表明,幂函数拟合效果最佳,且结合河网、河源密度与地形地貌确定阈值的方法是科学的;高翔等[6]基于30 m分辨率DEM,通过河网密度法确定祁连山国家公园流域河网阈值为6000个栅格时,提取结果相较于全国1∶250000三级水系数据中的河网更贴近真实水系。河网密度法可以解决阈值选取的主观性问题,但阈值大小仅影响模拟河网的疏密程度,而不会改变其空间位置。对于一些中小流域而言,DEM的分辨率较低(即栅格尺寸较大)时,断面宽度较小的河流所在的栅格可能会因同在一个栅格内的河流旁地形的平均作用,使得该DEM栅格的高程大于河流实际高程。提取河网时,此类DEM栅格可能会被识别为非河道栅格,从而造成提取到的模拟河网与实际水系位置存在偏差,故需结合区域实际影像对DEM进行预处理,然后利用河网密度法确定河源提取的最佳阈值,以增加河网提取的准确性。子流域的划分方面也开展了大量研究:Li等[7]基于高分辨率遥感影像,结合Canny算法提取流域特征信息,替代DEM数据构建河网,有效提升了太湖流域的子流域划分精度;Liu[8]等根据流域内的地貌特征,采用两级划分的方法对西藏内陆河流域进行划分,结果表明该方法在湖泊分类和子流域划分方面,能够有效地识别流域内的湖泊并保持各子流域的独立性;Stein[9]通过改进Pfafstetter系统,结合区域划分优化编码方法,在澳大利亚大陆实现了复杂河网的均匀子流域划分,显著提升了主干河道识别精度。众多学者就子流域划分的不同结果对流域径流等的影响也开展了大量研究[10-12]:Tripathi等[13]基于SWAT模型,分析认为不同子流域划分数量对印度东部南湾流域的土壤含水量等水平衡成分影响显著,但对地表径流的影响并不显著;Lin等[14]分析了子流域划分对西溪盆地的山坡泥沙生成模拟及其空间变化的影响,结果表明随着子流域划分阈值减小,山坡沉积物生成总量下降16.5%,空间变化增加97%;胡连伍等[15]研究了不同子流域划分结果对丰乐河流域径流、泥沙等的影响,提出不同的子流域划分层次对模拟结果的影响存在上下两个阈值,超过阈值则失真,低于阈值则难以反映真实空间分布。可见子流域的划分是分布式水文模型构建中十分重要的环节。目前,面向水文过程模拟的中小流域的子流域划分,主要是通过模拟的径流过程对子流域数量敏感性分析的统计学方法确定的[10-15]。当前研究在基于遥感影像定位实际水系,并结合数字高程模型DEM提取与之相符的流域河网方面已有一定进展。然而,利用此河网进行子流域划分后,关于子流域划分方案(数量与面积)对流域出口断面径流过程模拟精度影响的系统性研究仍显不足。
本文面向水文过程模拟,提出一种基于实际影像对主河道所在DEM栅格进行高程降低处理的方法,减少诸如圩区、洼地等地形地貌因素对于河网提取造成的偏差。采用河网密度法确定河网提取的最佳阈值[16],避免河网提取阈值选取的主观随意性。通过分析不同子流域集水面积阈值下的子流域划分结果,并结合我国一级小流域面积范围标准(50~300 km2),以及子流域出口均位于河网干流上的准则,确定子流域划分的合适阈值[17-18],比较分布式水文模型下不同的子流域划分方案在水文过程模拟中的差别。研究结果可为分布式水文模型的优化提供参考。
1 数据与方法
1.1 研究区域概况
丰乐河,古称桃溪,上接大别山区,下入巢湖,其源头分别为北支思古潭河、中支张店河和南支张母桥河,于龙咀合入丰乐河干流,下游汇入杭埠河后流入巢湖。桃溪流域(31°21′~31°47′N,116°26′~117°1′E)面积为1492.3 km2,属江淮丘陵区,为亚热带湿润性季风气候,多年平均降水量为1020~1200 mm,流域上游为浅山区,中下游多为丘陵、圩区,土壤主要为黄棕壤、潮土、水稻土等[19]。桃溪流域范围内高程最大值为468 m,最小值为-21 m,标准差为139.9 m,整体趋势为西部浅山区向东部丘陵区递减。
1.2 数据来源
DEM数据来源为地理空间数据云平台(www.gscloud.cn)的GDEMV2数字高程数据,分辨率为30 m。实际影像资料为天地图影像地图(www.tianditu.gov.cn)。水文过程模拟所用的降水、径流等资料由合肥水文水资源局提供。
1.3 DEM预处理方法
DEM的分辨率越高,越能真实地反映实际地形。然而,分辨率的提升会导致基于DEM的河网提取及子流域划分的工作量呈指数增加。同时,易出现很多集水范围较小的洼地[20],造成提取的河网栅格与实际水流方向存在偏差等情况[21]。通过降低实际水系主河道所在DEM栅格的高程,可使基于DEM提取的(等级较高的)河网水系与遥感影像中的水系更加一致。DEM预处理主要包括如下3步:
(1)根据实际水系与基于原始DEM提取的模拟河网,确定二者的偏差位置。在GIS中确定流域边界[22],然后根据实际遥感影像中的河道位置,提取流域边界内的实际水系(图1)。地理坐标系统一为WGS 1984。基于原始DEM提取出模拟河网,将模拟河网的主河道与提取到的实际水系进行叠加分析,识别主要的偏差位置。图2显示,自双河镇起下游区域模拟河网主河道与实际水系存在较大偏差。这表明,由原始DEM直接提取的模拟河网无法准确识别实际河流栅格。
(2)确定主河道所在DEM栅格的高程降低值。模拟河网主河道与实际水系存在偏差,可能是河流栅格受其附近地形的平均作用,导致河流栅格的平均高程值升高,河流栅格不能被准确地识别,而误将其周边的高程较低处识别为河道栅格;抑或河道周边存在圩区等较大面积的洼地,而洼地与河道之间的水流排泄建筑物未能在DEM中体现,易误将洼地等识别为河道栅格、实际河道识别为坡地栅格。故需依据影像提取实际水系的主河道,降低主河道所在DEM栅格的高程值,使流域内的产流都能通过主河道,进而通过出口断面流出。其中,降低河流所在DEM栅格的高程与填洼处理的功效类似,均为使流域内的产流能够通过出口断面流出,且降低河流栅格的高程还可促使提取的模拟河道与实际水系一致。经计算,实际水系主河道栅格在偏差处的平均高程为26 m,垂直于主河道方向上洼地栅格的最低高程为-13 m,主河道栅格平均高程与洼地栅格高程的最大差值为39 m。主河道栅格的平均高程降低后,应低于垂直主河道方向上洼地栅格的高程最低值,因此确定主河道处栅格高程降低值为40 m。DEM高程降低前、后的对比见图3。

图1桃溪流域实际影像提取水系
Fig.1Actual image extraction of Taoxi watershed drainage system diagram

图2原始DEM提取主河道与实际水系对比
Fig.2Comparison of extracted main stem and actual water system before DEM treatment

图3原始DEM(a)与高程降低后的DEM(b)
Fig.3Original DEM (a) and DEM after elevation reduction (b)
(3)降低主河道所在DEM栅格的高程值。进一步验证对桃溪流域实际水系主河道所在原始DEM栅格的高程值降低40 m后,提取模拟河网的合理性,其步骤如下:沿着实际水系主河道,以5 m作为步长增加栅格高程降低值,重复提取模拟河网并对比实际水系。当高程降低值达到40 m时,模拟河网主河道与实际水系完全吻合。进一步增加降低值,结果与40 m时一致,证明降低40 m为合理值。此时偏差处实际水系主河道栅格的平均高程为-14 m,低于垂直主河道方向上洼地栅格的高程最低值(-13 m)。将实际水系主河道所在DEM栅格的高程降低40 m,可迫使流域内的其他河流栅格流向实际主河道栅格,完成对原始DEM数据的预处理。DEM处理后提取的河流主河道与实际水系的对比见图4。
1.4 河网提取阈值确定方法
河网密度法通过设置一系列递增的集水面积阈值,生成对应的一系列模拟河网。计算不同阈值下的河网及河源密度,分别拟合河网密度、河源密度与集水面积阈值之间的关系曲线。随着阈值逐渐增大,等级较小的河流将被去除,当阈值增加到一定程度,将只保留主河道栅格。寻找关系曲线随集水面积阈值变化趋于平缓的点,即关系曲线拟合函数的拐点,该点对应的阈值为河源的最佳提取阈值。

图4DEM处理后提取的主河道与实际水系的对比
Fig.4Comparison of extracted main stem and actual water system after DEM treatment
1.5 子流域划分方法
采用不同的集水面积阈值划分子流域的具体步骤如下:
(1)先设置一个集水面积阈值(如10000栅格数)。(2)根据该阈值及由D8算法计算的流向和流量矩阵来确定流域内的河流网络。(3)沿着河流从流域下游向上游依次进行搜索,搜索出每段河流的起点,所有能流入该段河流起点的栅格编号记为该段河流的编号。如以流域出口作为起点的第1段河流编号为1,则全流域的栅格均能够流入到该段河流的起点,将全流域的河网栅格编号为1。再向上游搜索第2段河流起点,若该点所在栅格的上游集水面积超过设定的子流域集水面积阈值,则所有流入该起点的栅格编号改为2,以此搜索到上游每段河流,直至每个栅格都完成编号,栅格编号即为子流域编号。如此根据每段河流编号完成子流域编号,子流域内河流流向与其同编号的那段河流流向相同。(4)对于集水面积小于设定阈值的小面积子流域,则将其并入其流入的子流域,同时其编号也改为其流入的子流域编号,并且原小面积子流域后面的子流域编号依次减去1,完成修改。(5)将所有子流域编号完成后,就可以得出每个子流域范围,确定每个子流域的边界[23]。
本文确定最佳子流域提取阈值的方法,是通过设置一系列递增的子流域集水面积阈值生成子流域,计算不同阈值下的子流域数量及其对应的平均面积,得到子流域划分数量和平均面积与子流域阈值之间的关系曲线。随子流域集水面积阈值的逐渐增大,子流域划分数量有明显的减小趋势,子流域平均面积有明显的增大趋势;随阈值增加,当所有子流域的出口均位于干流上,且集水面积较小的子流域被合并时,此时的集水面积阈值即为最佳阈值。同时,尽量使划分的多个子流域面积的平均值处于我国一级小流域面积范围标准内,即50~300 km2。
1.6 水文过程模拟方法
采用三水源新安江模型对桃溪流域不同子流域划分方案下的日径流和次洪过程进行模拟。其中,子流域内的河网汇流采用滞后演算法,子流域出口至流域出口之间的河道汇流采用分段马斯京根连续演算法[24]。日径流过程模拟和次洪过程模拟的初始状态,均通过日尺度资料进行预热,预热期为1个月。日径流过程模拟以1982—2012年作为训练期,2013—2022年作为验证期,模型输入及输出资料均为日尺度;次洪过程模拟选取1982—2022年间20场最大洪水过程,以前15场作为训练样本,后5场作为验证样本,输入及输出资料的时间尺度均为6 h,以模拟流量过程与实测流量过程的确定性系数最大化为目标函数[25]。新安江模型所用参数个数为17个(包括流域蒸散发折算系数、上层张力水容量、下层张力水容量、深层蒸散发折算系数、流域平均张力水容量、张力水蓄水容量曲线指数、不透水面积占全流域面积比例、表层自由水蓄水容量、表层自由水蓄水容量曲线指数、表层自由水蓄水库对地下径流的日出流系数、表层自由水蓄水库对壤中流的日出流系数、地面径流消退系数、壤中流消退系数、地下径流消退系数、河网汇流滞时、马斯京根法河段传播时间和马斯京根法流量比重系数),各子流域之间产汇流参数取值不同,模型参数通过遗传算法率定。
2 结果与分析
2.1 河网提取阈值确定
由于河网提取过程中不同集水面积阈值提取的河网有一定的差别[26],尤其是河源的位置,且集水面积阈值的选取尚无明确的指导原则,具有一定的主观性。采用河网密度法,并以河源密度作为辅助验证,可有效避免阈值选取的主观性。分别以100~2400个栅格作为集水面积阈值,提取的流域面积为1492.3 km2,计算出不同阈值下的河网密度与河源密度,结果见表1。
表1桃溪流域水系特征信息
Tab.1Water system characteristic information of Taoxi Watershed

采用幂函数分别拟合河网密度、河源密度与集水面积阈值(栅格数)之间的关系,见图5。拟合函数变化趋于平缓的点即为最佳集水面积阈值点。

图5河网密度(a)、河源密度(b)与集水面积阈值拟合函数图
Fig.5Fitting function diagram of river network density (a) , river source density (b) and catchment area threshold
河网密度、河源密度与集水面积阈值(栅格数)的关系如下:

(1)

(2)
式中,N为河网密度(km/km2);n为集水面积阈值对应的栅格数;S为河源密度(个/km2);R2为确定性系数。
可见,两者拟合函数的确定性系数均接近1,拟合效果很好。再分别对式(1)、(2)求二阶导数(图6),寻找拟合函数的拐点对应的栅格数作为最佳阈值。当栅格数为900,即集水面积为0.81 km2时,河网密度与河源密度拟合函数对应二阶导数值几乎为0且后续数值几乎无变化,故选取该值为最佳集水面积阈值。
2.2 子流域划分
小流域作为水土保持治理的基础单元[27],可作为划分各子流域面积大小的重要依据。本文以我国一级小流域面积范围标准(50~300 km2)为子流域平均面积大小的参考,以每个子流域出口位于河网干流上为子流域划分的准则,划分桃溪流域的子流域。设置集水面积阈值从9 km2开始逐渐增加到144 km2(文中DEM栅格为30 m×30 m,即每10000个栅格对应9 km2集水面积,将集水面积阈值选为9 km2的倍数,便于确定子流域对应栅格数),计算出不同阈值下对应的子流域数量与平均面积。划分出的不同阈值下的子流域数量与平均面积见表2,子流域数量、平均面积与子流域集水面积阈值的关系见图7。

图6河网密度(a)和河源密度(b)拟合函数二阶导数图
Fig.6Graphs of the second derivative of the fitting function for river network density (a) and river source density (b)
表2不同阈值下划分子流域数量与平均面积
Tab.2The number and average area of sub-basins under different thresholds


图7划分子流域数量、平均面积与阈值的关系
Fig.7The relationship of the number of subbasins and their average area with the dividing threshold
将桃溪流域子流域集水面积阈值从9 km2逐步增加至108 km2。结果表明,划分的子流域数量由136个减少至7个,平均面积从11.0 km2增加至213.2 km2;当阈值从108 km2增加至144 km2时,子流域数量均为7个,平均面积均为213.2 km2。分析发现,当子流域划分的集水面积阈值较小时,划分出数量众多的小面积子流域(如阈值为9~36 km2时,子流域平均面积小于我国一级小流域面积范围标准),子流域的出口断面分布复杂,导致产汇流计算难度增加;而当阈值过大时,划分的子流域数量较少,难以充分反映不同区域的水文特性差异,从而影响模型对流域水文过程的真实模拟。进一步分析发现,当阈值为99和108 km2时,子流域数量分别为11和7个,平均面积分别为135.7和213.2 km2,均符合我国一级小流域面积范围标准,且子流域的出口均位于干流上。与划分11个子流域相比,阈值为108 km2时子流域数量减少了4个,说明流域内部分集水面积较小的子流域被进一步合并。在划分7个子流域的情况下,分布式水文模型不仅能够反映子流域内水文特性的空间异质性,还能显著减少计算工作量,提高模拟效率。对划分出的7个子流域进行空间拓扑关系编号(图8)。
分布式水文模型的构建与计算中,合适的子流域划分既有助于保持水文特性,又可优化计算效率。已有研究表明,子流域数量在一定范围内对径流模拟影响较小,而超过范围则会产生显著影响。例如,邱临静等[28]通过设置不同阈值划分杏子河流域子流域,采用SWAT模型模拟多年径流过程,该文将阈值从10 km2增加至100 km2的过程中,划分出子流域数量最低为7个,平均面积为216.4 km2,同时径流模拟结果没有较大变化;胡连伍等[15]应用SWAT模型,分析了丰乐河流域(该文提取的流域面积稍小于本文)子流域的不同划分结果对径流等模拟的影响,划分出子流域数量最低为5个,平均面积为257.8 km2。这些子流域的面积均在我国一级小流域面积范围标准(50~300 km2)内。本文中当子流域集水面积阈值为108 km2时,桃溪流域被划分为7个子流域,平均面积为213.2 km2,与上述文献中所划分的子流域的平均面积相近。尽管这些流域气候区和地形地貌存在差异,但划分出的子流域的平均面积却十分接近,表明子流域面积在我国一级小流域面积范围标准内的划分方法适用性较强,这一现象的潜在机理需要进一步研究。

图8子流域划分结果
Fig.8Division results of sub-watershed
2.3 日径流过程模拟
为验证上述对实际水系主河道所在DEM栅格的高程降低40 m后提取河网,以及将桃溪流域划分出不同数量的子流域对水文过程模拟结果的影响,设置3种子流域划分方案:方案A为将实际水系主河道所在DEM栅格的高程降低40 m后,按集水面积阈值108 km2划分的7个子流域;方案B为将实际水系主河道所在DEM栅格的高程降低40 m后,按集水面积阈值99 km2划分的11个子流域;方案C为由原始DEM,按集水面积阈值108 km2划分的7个子流域。基于1982—2022年的日降水、径流和蒸发资料,采用三水源新安江分布式模型[29]对3种方案下的日径流过程进行模拟,模型参数以实测流量和模拟流量的确定性系数最大为目标,采用遗传算法率定参数,其中每个子流域参数取值不同。以1982—2012年为训练期,2013—2022年为验证期。日尺度下的模拟结果为:方案A、B、C在训练期的确定性系数分别为0.87、0.87、0.86,在验证期的确定性系数分别为0.90、0.91、0.90。日径流过程模拟结果中,3种方案在训练期和验证期的确定性系数差别不大。包含最大洪水的2022年7月模拟和实测日径流过程见图9。

图9日径流过程模拟
Fig.9Simulation of daily runoff processs
2.4 次洪过程模拟
为进一步比较三者的模拟效果,基于桃溪流域时段长6 h的降水、径流等数据,选取1982—2022年间20场最大洪水过程进行次洪模拟分析。其中,前15场作为训练样本,后5场作为验证样本,次洪模拟结果见表3。训练期的15场洪水中,方案A、B、C的洪量相对误差在20%以内的均为14场,洪峰相对误差在20%以内的分别为15、15、14场,峰现时间误差在±6 h以内的分别为15、15、14场,平均确定性系数分别为0.93、0.94、0.93。验证的5场洪水中,方案A、B、C的洪量相对误差在20%以内的均为5场,洪峰相对误差在20%以内的均为5场,峰现时间误差在±6 h以内的分别为5、5、4场,平均确定性系数分别为0.94、0.93、0.94。根据国家标准《水文情报预报规范》(GB/T22482—2008)对20场洪水的合格率及精度等级进行评定。按照标准,每场洪水的洪量、洪峰预报相对误差绝对值小于许可误差(20%),峰现时间预报误差绝对值小于许可误差(6 h)为合格,方案A、B、C的合格场次分别为19、19、16场,合格率分别为95%、95%、80%。在20场洪水中,以确定性系数作为精度评定指标,方案A中甲级16场、乙级4场;方案B中甲级18场、乙级2场;方案C中甲级18场、乙级2场。选取4场洪水过程绘制径流过程线(图10)。

图104场洪水过程模拟
Fig.10Simulation of 4 flood processes
对桃溪流域3种子流域划分方案下的分布式水文模型进行计算分析,结果见表3。结果显示,对主河道所在DEM栅格高程进行降低处理的方案A和B,在20场次洪模拟中有19场合格,而未对DEM处理的方案C仅有16场合格。对于20200715号洪水,方案A、B的确定性系数分别为0.86、0.88,方案C的确定性系数为0.91;方案A、B的峰现时间误差均为-6 h,方案C的峰现时间误差为-12 h,方案A、B在峰现时间误差上优于方案C。结合降水资料分析可知,该场洪水期间桃溪站的降水强度大于上游雨量站,而流域面降水量通过泰森多边形法划分,方案C中桃溪站降水量在下游的子流域范围内占比相对方案A、B更大,因此,方案C中桃溪站降水显著影响下游汇流速度,导致洪峰提前。对于19840830号洪水,方案A、B的峰现时间误差均为6 h,方案C为18 h,由于流域上游雨量站的降水强度大于下游雨量站,且对主河道所在DEM栅格的高程降低处理后的方案A、B提取的主河道长度小于方案C,由此可见方案C从上游到下游汇流时间更长,因此计算的洪峰到达时间更晚。对于20100902号洪水,3种方案的洪量误差绝对值均超过20%。该场次洪水位于汛期末,且洪水前连续多日无降水,预热期计算出初始土壤含水量较低,导致整体洪量偏小,洪峰也均小于实测洪峰。综上,对主河道所在DEM栅格高程进行降低处理的方案A和B,对比未处理的方案C,在次洪模拟合格率上提升15%。同样对主河道所在栅格的DEM高程降低处理,划分7个子流域的方案A和划分11个子流域的方案B在20场次洪模拟的平均确定性系数分别为0.93和0.94。方案A虽然将更多小面积子流域进行合并,但在次洪模拟精度上并没有明显下降,说明划分为7个子流域是合理可行的。
表3桃溪流域20场洪水过程模拟结果
Tab.3Simulation results of 20 flood processes in Taoxi watershed

3 讨论
1)由于桃溪流域中下游丘陵区域的河道周边分布大量圩区,基于DEM栅格数据的河网提取中,高程较低的圩区栅格易被误识别为河网栅格。同时,由于河堤等建筑抬高主河道栅格高程,模拟河网与实际主河道存在偏差。文中高程降低处理方法在桃溪流域的适用性表明:对于地形更加平坦的平原地区来说,对实际水系主河道所在DEM栅格进行高程降低处理对河流栅格流向调整的影响更为明显,从而减少周边地形因素对河网提取的影响;而对于地形起伏较大的高山峡谷地区,由于不同河流栅格与周围的山地栅格高程数值相差较大,仅进行高程降低处理能否有效,仍需更多实例验证。同时,部分地区存在实际影像资料难以获取等问题,无法提取准确的实际水系位置,本文方法不再适用,如何增加本文方法在缺资料或无资料流域的适用性,值得进一步研究。
2)科学合理的子流域划分是流域分布式水文模型构建的基础,子流域面积的大小及其数量的确定尚无公认的准则。从水文过程模拟的角度看,无论是子流域所属的流域还是更大尺度流域的子流域,划分的最小子流域应保证该子流域范围内产生的地表径流、壤中流和地下径流均能从该子流域的出口断面流出,即子流域出口断面处的河床应切割到基岩,这样水文模型才能准确和方便地模拟出以子流域为单元、各子流域间通过河网建立水力联系的实际水文过程。目前,流域出口断面的河道深度以及河道下基岩埋深等相关数据难以获取。同时,本文方法主要面向水文过程模拟,仅从水文学角度进行计算,后续可以考虑结合行政区划边界,并增加从水力学角度考虑的因素进行子流域边界划分,该方面的研究值得拓展深入。
3)本文以我国一级小流域面积范围标准(50~300 km2)作为桃溪流域划分子流域平均面积的参考,并结合每个子流域出口均位于河网干流上的准则,与其他学者研究的一些子流域划分结果进行了对比,发现这些子流域的面积都在我国一级小流域面积范围标准内,但对于更多的中小流域而言,以我国一级小流域面积范围标准作为中小流域的子流域划分标准是否仍然适用,值得进一步研究。
4 结论
根据桃溪流域的实际影像资料提取实际水系,对原始DEM中主河道栅格进行适当高程降低处理,采用河网密度法寻找最佳集水面积阈值,再通过不同阈值划分子流域得到最合适的子流域数量,最后进行水文过程模拟比较,主要结论如下:
1)对主河道所在的DEM栅格进行高程降低40 m处理,可减少圩区、河堤等局部地形因素对河网提取结果的干扰。基于河网密度法,确定桃溪流域河源集水面积的最佳阈值为0.81 km2。通过数理方法确定集水面积阈值,有助于避免阈值选取的主观性,提取出的模拟河网可有效减少伪河道的生成。
2)以108 km2作为子流域集水面积阈值,划分出7个子流域,每个子流域的出口均位于河网干流上。子流域平均面积为213.2 km2,划分结果与我国一级小流域面积范围标准(50~300 km2)一致,也与一些其他地区划分的子流域平均面积接近。
3)分别进行DEM高程降低后划分7个子流域、高程降低后划分11个子流域和原始DEM划分7个子流域3种不同方案下的新安江模型水文过程模拟,日尺度下3种方案在训练期确定性系数分别为0.87、0.87、0.86,验证期分别为0.90、0.91、0.90,三者精度差别不大;20场洪水过程模拟中,3种方案的合格率分别为95%、95%、80%,前两种方案的合格率明显高于第3种方案,而前两种方案次洪模拟的平均确定性系数分别为0.93和0.94,差别较小。因此,将桃溪流域划分为7个子流域有助于提高次洪模拟合格率,减少分布式水文模型计算的工作量。