郑佳玮，陈兴伟，2,3*，李孝成，阮伟芳

(1.福建师范大学地理科学学院，福建 福州 350007；2.福建省陆地灾害监测评估工程技术研究中心，福建 福州 350007；3.湿润亚热带山地生态国家重点实验室培育基地，福建 福州 350007；4.福建省水利水电科学研究院，福建 福州 350001)

流域内所有河流、湖泊等各种水体组成的水网系统，称为自然界的水系，其中河流是水系最为重要的组成部分[1]。水系的提取通常是根据局部地形特征，借助一定算法，辅助获取研究区域内的水文要素信息[2]。数字高程模型(Digital Elevation Model，DEM)是一组有序数组，用以描述地面高程值空间分布的实体地面模型，能够在一定空间分辨率下反映地区局部的地形特征。因此，成为目前用于水系提取分析的主要数据[3]。基于DEM数据提取的水系是对真实水系的概化表达，最大程度上反映了水系的真实特征。应用DEM提取水系主要包括2个步骤，一是流向计算[4-6]，目前根据水流方向计算的模型可以分为单流向算法(SFD)和多流向算法[7](MFD)。单流向算法以D8算法最为典型[8]，在大尺度区域河网提取中，D8算法应用最为广泛[9-10]。

应用DEM提取水系的第2个步骤是水系精细程度的确定[11]，除了王林等[12]采用最小河流长度以确定提取水系精细程度的方法外，目前主要采用控制集水区域面积的大小来确定提取水系的精细程度。集水区域面积越小，提取的水系越密集；反之，提取的水系就越稀疏。因此，理论上存在一个集水面积阈值。集水面积阈值一般定义为支撑一条河道永久性存在所需的最小集水区域面积[13-14]。Giannoni等[15]提出根据坡地水流路径的频率分布法确定集水面积阈值，孔凡哲等[16]提出根据集水区域面积与河网密度关系曲线寻求最佳阈值；李丽[17]将实际水系的分形维数与不同集水区域面积下的水系分形维数进行比较，得出阈值；吴泰兵等[18]提出基于改进适度指数法确定集水面积阈值；Jefferson等[19]总结出形为AxSy的流域内坡度、集水面积尺度化指数，即给定阈值，超过该阈值的栅格属于河道；张鹏举等[20]通过河道分支比确定集水面积阈值。上述方法普遍利用集水区域面积与水系特征量的变化曲线来确定集水面积阈值，通过目视判断水系的集水面积阈值，存在一定的主观性。杨邦等[21]采用1∶82万纸质地图的水系；程中阳等[22]以1∶25万的全国水系图作为实际水系，选用河道平均坡降法、河网密度法、流域宽度分布法、分形维数法等方法计算集水面积阈值，讨论各方法的适用性。这种提取水系与实际水系对比方法的主要问题是何种比例尺水系适合作为判别标准及其与DEM数据分辨率的匹配问题，如何定量比较实际水系和提取水系之间差异大小的方法，也仍需进一步研究。

陈冬平等[23]提出利用提取的水系河源与实际河源的距离误差最小的方法确定集水面积阈值，以避免目视判断集水面积阈值的主观性。河源即河流发源地,指该河干流离河口处最远开始补给的地方[24]。理论上，提取水系的河源应该与实测河源重合，而集水面积阈值的合理确定，可以准确反映提取河源的位置，提取水系的河源能否合理反映实测河源的位置，可以作为一种判断提取水系合理性的评价标准。

因此，本研究以河源的合理提取为评判标准，以东南沿海的晋江流域为研究区，选用基于DEM提取水系的河网密度法、分形维数法、河源最小误差法和河道分支比法等4种主要方法，以比较分析各方法提取水系的差异，探讨各方法提取水系的合理性。

1 数据来源与研究方法

1.1 研究区概况

位于福建省东南部的晋江，是该省的五大河流之一。晋江流域所处的地理范围在117°44′～118°47′E，24°31′～25°32′N，位于亚热带海洋季风气候区，多年平均降水量为1 147 mm。流域气候暖热，冬短夏长，热量丰富，日照充足，热量资源丰富。降水较多而时空分配不均，春夏多雨，秋冬少雨，降水量的年际变化也较大。台风、暴雨、低温、干旱等灾害性天气时有发生。晋江发源于福建省中部戴云山，河源位置位于泉州市安溪县境内桃舟乡达新村梯仔岭东南坡[25]，坐标位于117°41′42″E，25°21′30″N，河流全长302 km，总面积5 269 km2。晋江上游分为东溪和西溪两大支流，沿途流经永春、安溪、南安、晋江、鲤城、丰泽等县市区(图1)。

图1 晋江流域位置

1.2 数据来源

研究区数据是由美国航空航天局(NASA)和国防部国家测绘局(NIMA)联合测量的STRM DEM数据(http://dwtkns.com/srtm30m/?tdsourcetag=s_pcqq_aiomsg)，分辨率为30 m，数据覆盖全球陆地表面的80%以上，该测量数据覆盖中国全境。

采用ArcGIS 10.2 Hydrology模块的相关功能，通过以下步骤提取水系：①运用Fill工具进行DEM洼地填充；②采用基于D8算法的水流方向计算工具，利用Flow Direction工具完成水流方向的计算；③运用Flow Accumulation工具计算集水区域面积；④通过设定不同的集水区域面积，得到不同密度的水系图。

1.3 研究方法

1.3.1河网密度法

河网密度是指流域内干支流的总河长∑L(km)与流域面积A(km2)的比值，基本单位为(km/km2)。河网密度的计算公式为：

(1)

通过设定不同的集水区域面积，生成一系列模拟的河网，计算各模拟河网的河网密度，得到模拟河网密度随集水区域面积的变化关系曲线。河网密度变化趋于平缓时对应的集水区域面积即为集水面积阈值[16]。

1.3.2分形维数法

目前有2种方法广泛用于水系分形维数的估算，一是基于分形定义的盒维数法，二是基于Horton定理。其中，基于盒维数法对水系的分形维数估算可以适应不同类型的线性体，有助于分析计算复杂河网的形态[26-27]。

盒维数法以边长为r的栅格覆盖分形体，将与分形体相交的格网计为非空，否则给出空值，并统计非空栅格数计为N(r)。当r持续变化时，即对应一组r1，r2，r3，……，得到一组N(r1)，N(r2)，N(r3)……。以点[lgr，lgN(r)]为坐标作散点图，若测度r与非空格网数N(r)在双对数坐标轴上呈现线性关系，那么用最小二乘法可以拟合一条直线:

lgN(r)=-DHlgr+C

(2)

式中DH——水系的分维值；r——栅格边长；N(r)——包含河流片段的非空栅格数；C——待定常数[28-30]。

绘制不同集水区域面积水系的分维值随集水区域面积的变化关系图，水系分维值变化趋于平缓的点所对应的集水区域面积即为集水面积阈值[11]。

1.3.3河源最小误差法

河源最小误差法采用有关部门综合考察确定的河源来评估提取水系的合理性。受地形、算法等其他因素的影响，设定不同集水区域面积提取的河源与实际河源会存在一定的距离误差。运用实际河源与提取水系河源两点之间的距离误差作为判断水系提取合理性的标准，将河源距离误差最小时的集水区域面积看作是集水面积阈值[23]。

1.3.4河道分支比法

理论上，真实水系河道分支比的均值趋于一个常数。虽然河道数目随集水区域面积的变化而发生变化，但各级河道数量比几乎保持不变，河道分支比的平均值在某一范围内趋于稳定。据此，可以绘制平均河道分支比随集水区域面积的变化关系图，以河道分支比变化趋于稳定时的集水区域面积，作为集水面积阈值[20]。

1.3.5集水面积阈值判定的一阶导数法

河网密度法、分形维数法以及河道分支比法都是基于计算的水系特征量，绘制与集水区域面积的变化关系曲线，认为关系曲线趋于平缓时所对应的点即为集水面积阈值。由于难以用目视的方法判别关系曲线趋于平缓的点，可用关系曲线的一阶导数判定集水面积阈值。分别对上述水系提取方法得到的河网密度、分形维数、河道分支比与集水区域面积的变化关系曲线进行曲线拟合，进而求得拟合曲线的一阶导数，判定集水面积的阈值。

2 结果与分析

2.1 河网密度法

依次设定不同的集水区域面积，求得河网密度与集水区域面积的变化关系见图2。当集水区域面积(即栅格数)增加，河网密度呈现下降的趋势，下降的速率由急剧变为缓慢，当栅格数从0增加到10 000时，河网密度下降较为剧烈，之后趋于平缓。对河网密度拟合函数求一阶导数，在10 000时一阶导数趋于平缓且接近0，此时的集水区域面积为9 km2，即为河网密度法确定的晋江流域集水面积阈值。

2.2 分形维数法

集水区域面积与水系分维的关系变化及其拟合函数曲线的一阶导数见图3。结果表明，当栅格数达到18 000左右时，一阶导数曲线趋于平缓且接近0，则分形维数法确定的晋江流域集水面积阈值为16.2 km2。

图3 分形维数与集水区域面积变化关系

2.3 河源最小误差法

在500～24 000栅格数范围内，选取系列栅格数生成不同密度的水系图。在ArcGIS 10.2下加载晋江河源数据与提取的水系图，使用measure工具量测提取水系河源与实际河源坐标的距离误差，从而得到距离误差与栅格数的变化见图4。从图4看出，当栅格数增加到6 000个左右时，距离误差由大到小；而当栅格数从6 000个左右继续增加时，距离误差又呈现由小到大的变化。为了进一步确定最小误差及其对应的栅格数，通过函数拟合，得到曲线的方程为：

图4 河源距离误差与集水区域面积变化关系

y=1E-05x2-0.173x+1526.1

(3)

上式一阶导数等于零时的栅格数为6 315，此时提取河源与实际的距离误差最小。相应地，河源最小误差法确定的晋江流域集水面积阈值为5.68 km2。

2.4 河道分支比法

通过设定不同的集水区域面积得到不同密度下的晋江流域河道分支比，所得河道分支比与集水区域面积的变化关系及其拟合函数曲线的一阶导数见图5。图5表明，当栅格数达到12 000时，河道分支比拟合函数的一阶导数趋于平缓且接近0，对应的10.8 km2被认为是提取水系的集水面积阈值。

图5 河道分支比与集水区域面积变化关系

3 讨论

3.1 集水面积阈值确定的不确定性

河网密度法、分形维数法以及河道分支比法皆是基于河流形态变化趋于稳定来确定阈值，对阈值的判定主观性较强，存在较大的不确定性。而河源最小误差法通过实测河源与实际的距离误差和集水区域面积的变化关系，拟合出的函数曲线是一元二次函数，具有先减小后增大的特点，存在最小值，即用距离误差最小时所对应的集水区域面积作为阈值，保证了阈值确定的合理性。

3.2 水系提取结果合理性分析

4种集水面积阈值确定方法得到的晋江流域水系分布见图6。4种方法提取的水系支流的形态特征变化较为明显，集水区域面积影响提取水系信息的详细程度。基于分形维数法和河道分支比法确定的集水面积阈值较大，提取出的水系形态结构较为简单，而基于河源最小误差法以及河网密度法确定的集水面积阈值较小，提取出的水系结构较为复杂，河道数目较多，结构较为精细。

理论上，4种方法提取河源与真实河源应该重合，实际上，提取河源与真实河源之间存在距离，称为距离误差。采用ArcGIS 10.2的Measure测量工具对提取河源与实测河源进行距离测量，结果为距离误差。4种方法提取的晋江流域河源的距离误差见表1。由于算法、数据精度等因素的影响，4种方法提取的晋江流域河源与实际河源都有一定误差，但河源最小误差法的误差最小，所提取的水系应能更合理反映晋江流域的真实水系。

a)河网密度法

c)河源最小误差法

表1 不同方法的集水面积阈值及其距离误差

从原理上看，河网密度法、分形维数法和河道分支比法是从河流形态的角度探究水系的形状结构特征，以找到河流形态特征趋于稳定的集水区域面积，但由于判断稳定与否的指标不尽合理，阈值的确定方法有待进一步改进。而河源最小误差法引入实际河源位置来控制水系的提取，一定程度上保证了提取水系的合理性。实际上，河源最小误差法将提取水系与实际水系的对比，简化为提取河源与真实河源的对比。河源即是水系起点，由于水系自身形态受地形控制，合理确定了水系的起点，水系本身误差应该也比较小。此操作方法简便，在难以获取实际水系数据的情况下具有较高的可行性。

4 结论

以晋江流域为例，选用河网密度法、分形维数法、河源最小误差法和河道分支比法等4种方法，分别确定水系提取的集水面积阈值，提取晋江流域水系；以实际河源为基准，比较不同方法提取河源与实际河源的距离误差，分析4种水系提取方法的合理性，结果表明。

a)水系提取的河网密度法、分形维数法和河道分支比法，其集水面积阈值确定的目视判别，具有较大的不确定性；这些方法提取的河源，与实际河源的距离误差都较大。

b)河源最小误差法提取的河源与实际河源的距离误差具有驻点，据此确定集水面积阈值所得到的距离误差最小，以此确定的河流水系具有唯一性，从而保证了提取水系的合理性。因此，河源最小误差法是晋江流域4种基于DEM真实水系提取方法中比较合理的方法。

本文工作仅以东南沿海的晋江流域为例开展研究，所得结论的推广还需进一步考虑不同地貌特征流域的适用性、不同类型和分辨率DEM数据可能的影响等问题。