1 引 言

汇水区域是指地表径流或其他物质汇聚到一共同出水口的过程中所流经的地表区域^[1]，是河流级别、长度、分布密度等几何特征的综合反映^[2]。在地图水系综合中，汇水面积是河流选取的关键量化指标^{[3, 4]}。目前基于河流数据的汇水区域提取均采用“空间均衡竞争”的策略，主要有两种途径：① 栅格提取方法，该方法将各线状河流栅格化，通过同步膨胀运算获取河流的Voronoi图，可近似得到每条河流的汇水区域范围^{[4, 5]}；② 矢量提取方法，该方法利用河网矢量数据构建Delaunay三角网，提取相邻河段间的三角网骨架线，从而建立每条河段的汇水区域多边形并最终求得河流汇水区域^[2]。上述方法都是从河流单一要素出发，没有考虑地形因素的作用。由于现实地形复杂多变，譬如两条河流之间为强烈非对称地形时，按照现有方法提取的汇水区域可能与实际情况存在较大偏差。如图 1所示，实线a为实际山脊分水线，虚线b为按“空间均衡竞争”方法划分的河段1与河段2之间的汇水分界线，两者之间存在较大差异，处理不当将影响到河流选取的结果。

图 1 非对称地形下的河流汇水区域示意图 Fig. 1 The asymmetrical mountain between river segments

图选项

另一方面，单独基于地形信息的汇水区域提取多从规则格网DEM数据出发，采用D8、DEMON、Dinf等路径算法^{[1, 6]}。此类方法将提取的所有谷底线均视为河流，只考虑地形信息而未考虑真实的河流状况，往往导致两个问题：① 河流是特定的地理目标，沟谷并非都是常年有流水的河流；② 由于局部平坦地形的存在，实际河流与提取的谷底线并不完全一致。因而此类方法不适用于解决水系综合中的河流选取问题。

由于河流与等高线具有天然的空间耦合关系，即河流常被看成是地形的“骨架”，而地形构筑河流的流域范围^[5]，因此要更准确地反映汇水区域需要同时考虑河流和地形两个要素的特征信息^{[1, 6]}。本文结合矢量河网与等高线数据，提出了基于双要素协同的平三角形骨架线提取规则与梯度向量引导的分水线搜索规则，然后联合使用两类规则提取分水线并建立各河段的汇水区域。

汇水区域提取的关键在于汇水边界线即分水线的提取。由于山脊线是天然的分水线，因而汇水边界线可主要通过山脊线提取获得，但由于山脊线并非连续且封闭的曲线，仅通过山脊线的提取无法获得完整的河流汇水范围^{[7, 8, 9, 10, 11]}。因此，本文提出了两类分水线：显式分水线与补齐分水线。提取这两类分水线即可构建封闭的汇水区域多边形。

显式分水线是从等高线山脊弯曲段上按一定规则直接提取的山脊分水线(如图 2中实线所示)，是山脊点的有序连线。

图 2 显式分水线段与补齐分水线段的提取 Fig. 2 Extraction of explicit and supplemental watershed segment

图选项

补齐分水线是为了使得汇水边界线连贯、完整，而按一定规则连接显式分水线与汇水点或地图边界所形成的线段，如图 2中虚线所示，其中O点为汇水点。

汇水点是支流向主干河流汇水的交点，从各个汇水点出发沿相邻河段之间的山脊区域提取分水线，可以解决传统方法提取的分水线破碎且连接困难的问题，且能避免当前尺度下细微地形起伏对分水线提取的噪音影响。

等高线与河流曲线统称为约束线，对于任意三角形边，若边的两个端点位于同一条约束线上则称这两点共约束线，若边在约束线上则称此边为约束边。

分水线段是指从约束Delaunay三角网(CD-TIN)三角形中提取出来的分水线线段或分水线线段有序连接形成的弧段。

本文通过构建CD-TIN来对等高线与河流双要素数据进行融合，三角形特征是判断成组三角形组合模式的基础，也是进行分水线段提取的基本条件。根据顶点所在约束线要素类型与顶点间邻接关系的不同，可将CD-TIN中的三角形分类，表 1为与分水线提取相关的三角形分类及其形式化描述。

CD-TIN中的顶点集合可表示为V={v₁,v₂,…,v_k}，三角形边可表示为e={<v_a,v_b> | v_a,v_b∈V }。对于任意给定的三角形t(v₁,v₂,v₃)(简称为Δt，其中v₁、v₂、v₃为3个顶点)，定义CLNum为Δt中共约束线的顶点对数；SeqNum为Δt中的约束边个数；JointNum表示Δt中汇水点个数，一般为0或者1；定义TypeNum(ftype)为Δt中包含所在约束线要素类型为ftype的顶点个数，参数ftype∈{Contour,River}，Contour、River分别表示等高线与河流类型。各类三角形分布位置如图 3所示。

图 3 CD-TIN三角形分类示意图 Fig. 3 Classification of CD-TIN triangles

图选项

将双要素自身及其之间叠加形成的顶点集合采用逐点插入算法构建Delaunay TIN，采用文献^[12]的算法嵌入河流与等高线约束线建立CD-TIN，最后建立其拓扑关系，并存储相关的地形信息与属性信息。

实际地形中的河网已大致将整个研究区域分成了多个未封闭的子区域，若从各个汇水点出发，基于CD-TIN特征三角形迭代提取并连接显式分水线与补齐分水线，依此对子区域进行细分，可最终形成包围各河段的封闭的汇水区域。

显式分水线构成了汇水边界的主要部分，本文提出了基于骨架线^[13]的提取规则与梯度向量^{[7, 14]}引导的分水线搜索规则。骨架线提取规则分别按特定的规则提取不同类型三角形内部的骨架线，并将邻接的三角形骨架线连接形成连续的分水线段。梯度向量引导的分水线搜索在三角形内完全或部分沿着搜索起点的梯度方向增长分水线。显式分水线提取需要根据当前所在三角形特征选择合适的提取规则，在有平三角形的等高线山脊弯曲部位可按骨架线规则提取，如图 4(a)；在无平三角形的等高线山脊弯曲部位则按梯度算法进行提取，如图 4(b)。

图 4 不同地形条件下的分水线提取示意图 Fig. 4 Extraction of watershed in various terrain conditions

图选项

补齐分水线的提取也分两种情况：一是连接汇水点与邻近等高线山脊点的路径提取，如图 4(c)虚线所示(部分)，可根据汇水点邻接的三角形组合类型按规则进行提取；二是连接显式分水线与区域边界的路径提取，如图 4(d)虚线所示，可按规则适当延长显式分水线到达地图边界。下文将详细描述各种分水线提取规则。

骨架线规则适用于山脊区平三角形与河流汇流区相关三角形的特征线提取。山脊部位构建CD-TIN后，在等高线弯曲凹部一般会形成平三角形簇，而且在两相邻等高线弯曲之间会形成桥接三角形。对具有此类特征的区域，可根据各类三角形特征分别提取骨架线作为分水线段，具体提取规则如下。

规则1：侧翼平三角形取平行约束边的三角形中位线作为分水线段，如图 5(a) ΔA。

图 5 等高线弯曲部位平三角形特征线提取 Fig. 5 Extraction of feature line in the plain triangles at the contour bend

图选项

规则2：突部平三角形取过突顶点(等高线上两约束边的公共顶点)的三角形中线作为分水线段，如图 5(a) ΔB。

规则3：桥接三角形取两顶点共约束线的边的中线作为分水线段，如图 5(a) ΔC。

规则4：为避免分水线出现较大的锯齿波动，将侧翼平三角形内的骨架线连接成分水线段后进行化简处理，使得分水线平直。

规则5：分岔平三角形可分为3种子类型，设最长边与最短边之比为r，若r≥2则判定其为Ⅰ型；若r≥1.5且r ＜2则判定为Ⅱ型；r ＜1.5则判定为Ⅲ型。Ⅰ型取平行于最短边的中位线作为分水线段，Ⅱ型取最长边上的中线作为分水线段，Ⅲ型取重心到各边中点的连线作为分水线段，如图 5(b)。

受等高线顶点密度的影响，并非所有等高线弯曲部位都可形成平三角形，按照骨架线规则无法提取到所有的分水线段。依据流水总是沿斜坡最陡方向流动的原理，本文提出了梯度向量引导的分水线提取算法(简称梯度算法)，由低至高追溯获得分水线段，比较适用于无平三角形斜坡上特征线的提取，且能使得连接分水线段后的分水线保持连贯、完整。

设分水线搜索起点为B，梯度向量(用G表示)按文献^[7]方法进行计算，以B为起点的梯度方向射线与三角形边或边延长线的交点记为C_i，梯度算法的提取流程可描述如下：

(1) 对于搜索起点B，若为V(CD-TIN顶点集合)中的顶点则获取其高程值h，然后查找与点B连接且高程大于h的三角形(其高程取重心高程值)，存入集合T₀并转步骤(2)；否则，若点B非V中的顶点则查找与点B所在边邻接的较高处三角形，然后转步骤(3)。

(2) 若T₀中存在邻接平三角形的三角形，则将其作为待处理三角形并转步骤(3)；否则，取T₀中高程值最大的三角形子集T₁，若T₁中元素唯一则转步骤(3)，反之若有多个元素则分别计算梯度向量并求取交点C_i，取T₁中交点C_i在内部的三角形为待处理三角形并转步骤(3)；若以上条件皆不满足则获取T₀中各三角形上与点B邻接的顶点集合V₀，然后取V₀中离点B最近且高程大于h的顶点增长分水线并转步骤(4)，如图 6(a)。

图 6 沿梯度向量搜索分水线(箭头为高程升高方向) Fig. 6 Extraction of watershed directed by gradient vector (arrow points to the higher place)

图选项

(3) 对于待处理的高处三角形t_h，若与平三角形邻接则直接连接平三角形顶点并转步骤(4)，如图 6(b)；否则计算交点C_i，并连接点B与点C_i增长分水线，如图 6(c)，然后将搜索起点B更新为交点C_i并转步骤(4)。

(4) 判断是否满足终止条件：① 与其他分水线会合；② 到达平三角形或地图边界；③ 超过生长步数阈值(可能出现增长异常)。满足条件则终止提取过程，否则转步骤(5)。

(5) 迭代执行步骤(1)至步骤(4)，逐阶向较高等高线逼近并增长分水线，直到终止。

由于河段上的顶点不具备高程信息，无法采用显式分水线提取方法来获得补齐分水线。而与汇水点邻接的多个三角形往往具有典型的组合形态，汇水点邻接三角形以河段为边界可形成3到多个群组，群组组合模式有两类：带状模式与扇形模式。带状模式是指群组内的三角形互相邻接而成条带状组合分布，如图 7中河段R₀与R₁之间的三角形群组；扇形模式是指群组内的三角形互相邻接而呈发散的扇形分布，如图 7中R₁与R₂之间的三角形群组。针对以上两种组合模式，可按规则6提取其内部分水线段。

图 7 汇流区三角形内分水线提取 Fig. 7 Extraction of supplemental watershed in the triangles near river confluence regions

图选项

规则6：对于带状模式，根据各三角形类型分别提取内部特征线，如图 7中虚线Ov₅所示，提取规则与平三角形骨架线提取类似。对于扇形模式，找出三角形群组中在河流约束线上的约束边，在图 7中即为Ov₃、Ov₆；计算两边之间的夹角并求其角平分线，然后查找边集合中与角平分线所构成夹角最小者作为分水线段。

规则7：当提取显式分水线到达边界三角形时，可将上一段分水线段向地图边界方向延长直至与地图边界相交，取延长线段作为补齐分水线段。

由河流数据计算获得所有河流汇水点，分别从每个汇水点出发提取分水线段。对于给定的汇水点O，交于O点的河段一般为3条(在复杂地形中可能多于3条)，理论上也应有3条分水线交于点O。从点O出发提取分水线的算法过程可描述如下：

(1) 搜索点O的一阶邻接三角形集合T₀，并以点O为中心按顺时针进行排序，然后以交于点O的河段为边界将T₀划分为子集{T_a,T_b,T_c,…}，从每个子集所在的分区出发分别提取一条分水线。

(2) 从点O出发，按规则6提取补齐分水线段，到达分区内离点O最近的等高线顶点A(如图 2线段OA)。

(3) 以A为起点，按梯度向量搜索规则向高处增长分水线，直到有平三角形存在(如图 2点B)，则转到下一步。

(4) 以B为起点，按骨架线提取规则提取平三角形(包括桥接三角形)内的显式分水线段，然后向较高位置的平三角递进增长，至无法搜索到邻近平三角形时再转到步骤(3)按梯度方法增长。

(5) 根据条件交替执行步骤(3)、(4)直至到达地图边界区域(如图 2点C)或检测到其他已提取的分水线时(如图 2点E)终止提取。

(6) 根据提取终止情况，按原分水线增长方向延长，并连接其与边界线或已有分水线的交点，即可获得由补齐分水线段与显式分水线段组合构成的一条连续的分水线。如图 2，弧段OA、CD为提取的补齐分水线段，弧段AB、BC为显式分水线段。

经过处理，整个研究区的区域边界、分水线段及其交点共同构成了分水线网络(如图 8)。其中网络节点集由汇水点集、分水线连接点集、分水线与区域边界交点集等3个子集组成，连接网络节点间的分水线段即可获得分水线网络的边集合。

图 8 分水线网络示意图 Fig. 8 Diagram of the watershed network

图选项

分水线网络可将分水线数据进行结构化组织并存储其拓扑关系，结合已有的河段拓扑关系，可采用多边形拓扑构建算法自动生成各汇水区域多边形。根据已有的CD-TIN拓扑信息和顶点的空间索引，可快速建立CD-TIN中各个三角形与汇水区域之间的关联，对汇水区域内三角形的面积求和即可得到所属河流的汇水面积。

本文采用陕西省绥德地区1∶5万DLG数据，以Visual C#开发程序进行了汇水区域提取试验。源数据为一组单线河流与等高线，如图 9(a)，其中等高线28条，河流27条，总点数2906个，构建CD-TIN后三角形为5776个，如图 9(b)。按本文方法，提取的封闭汇水多边形，如图 9(c)，红色曲线为河段的汇水区域边界线。

图 9 河流汇水区域多边形的构建过程 Fig. 9 Construction of catchment basins of the rivers

图选项

本文采用双要素协同的提取算法、河网Delaunay三角网骨架线算法与栅格Voronoi图算法(简称方法1、方法2与方法3)分别进行了汇水区域的提取试验，图 10(a)中的红色实线为方法1提取的汇水区域边界线，其中的数字为河段序号，图 10(b)、(c)分别为方法2、方法3的提取结果，图 10(d)为3种方法提取结果叠加效果图。通过对比，后两种方法提取结果比较相似，与等高线山脊弯曲部位吻合较差，局部区域明显偏离山脊而靠近河流；而本文方法提取的分水线与地形山脊吻合较好。

图 10 3种算法提取结果对比 Fig. 10 Comparison of catchment basins extracted by three algorithms

图选项

表 2为3种算法提取的汇水区域面积结果，其中ID表示河段序号，S_A、M_A、V_A分别表示方法1、方法2与方法3的汇水区域面积计算结果，γ_SM表示M_A与S_A的相对误差值，γ_SV表示V_A与S_A的相对误差值，S_O、M_O、V_O分别表示3类算法提取结果按面积降序排列的序号。通过分析可以发现，后两种方法提取结果与方法1提取结果相比，平均相对误差值高于25%。河段汇水面积越大表示该河段重要性程度越高，第15号河段通过3种算法计算排序结果分别为5、12与10，通过方法1计算，该河段在河网综合选取中的重要性得到较大提升，将会被优先保留而使得综合结果更为合理。

本文基于CD-TIN来统一组织双要素的地形信息并对三角形进行分类，然后根据各微观地形特征分别按骨架线算法与梯度算法提取了显式分水线与补齐分水线，对各分水线段进行连接形成分水线网络，并在此基础上按拓扑多边形构建算法实现了各河段汇水区域多边形的构建。较之单纯基于河网的三角网骨架线方法与栅格Voronoi图的提取方法，本方法更多地考虑了地形信息在流域提取中的作用；与基于DEM的提取方法相比，本方法从矢量数据中提取汇水区域，提取结果可直接使用，并避免了基于栅格数据的提取结果在矢量化后存在的边界不明确问题。本文方法适用于水系与等高线具有较强协同关系的各种地貌类型下的汇水区域提取。下一步的研究将结合水系树建立多尺度的流域层次结构，并基于此进行地貌与水系协同的自动制图综合。