地貌与水系的自动综合是地图自动综合的重要研究内容，而水系作为地貌的骨架线，它的综合对其他要素有一定的控制和制约作用，其重要性通常要高于地貌要素^[1]。近年来，关于水系的综合研究众多，主要集中在水系选取方法上，取得了较多的成果^{[2, 3, 4]}，而对水系曲线化简的研究较少。究其原因，是因为水系各河流之间的等级结构使其差异明显，更适合于选取，并且通常也认为水系的形态相对简单，采用常规的曲线化简方法即可。

目前，曲线化简方法分为以下3类：一是基于几何特征分析的曲线化简方法，这类方法也称作点位信息度量方法^[5]，即从几何特征角度分析曲线上各点的信息量大小，从而决定点是否选取，这类方法主要有D-P法^[6]及其改进方法^{[7, 8, 9]}、LI-OPENSHOW法^{[10, 11]}；二是基于结构特征分析的曲线化简方法，这类方法以曲线的弯曲结构为分析对象，通过对弯曲进行提取、结构组织与特征分析从而确定弯曲的保留与删除^{[12, 13, 14, 15]}；三是基于数学分析的曲线化简方法，主要有采用分形分析方法^{[16, 17]}和小波分析方法^{[18, 19]}。

然而，以上曲线化简方法都是针对二维曲线，且是对单条曲线处理。因此，当应用到水系曲线化简时，就存在以下两个问题：一是化简难以体现三维形态特征，特别是面对水系的“叠水”现象，难以保留其“叠水”特征(如图 1所示)，虽然文献^[20]提出利用3D D-P算法对等高线进行三维综合化简，但该方法的本质是将等高线看成三维散点，进行“体”的综合，并不是直接对三维曲线进行化简，故其并不适合于水系三维曲线化简；二是难以维护河流之间的拓扑关系，常规处理方法是将河流数据在交汇点处先分段再化简，如文献^[21]采用分段D-P化简方法保留交汇点从而维护其拓扑关系，然而，河流分段将破坏河流曲线的整体形态特征。事实上，要保持水系拓扑关系的一致性，并不意味着水系交汇点的几何位置在化简过程不能变化，反而，从多尺度表达的角度上看，水系交汇点的几何位置原本就应该随着水系综合的程度而相应变化。

图 1 水系的“叠水”特征 Fig. 1 The superimposed characteristic of water system

图选项

因此，本文对D-P算法进行三维曲线化简扩展，提出一种顾及拓扑一致性的水系三维曲线化简方法，在化简过程中既能保持水系三维形态，又能保证水系拓扑关系。

顾及拓扑一致性的水系三维曲线化简方法，需要解决两个问题：一是如何保持三维形态特征；二是如何维护水系中河流之间的拓扑关系。针对前者，将D-P算法进行三维曲线化简扩展，以点到空间直线的三维点线距代替平面上的点线距，从而使特征点的选取能顾及到整个三维形态特征。对于后者，通过构建水系树结构表达河流拓扑关系，并按照先干流后支流的顺序进行化简，在化简之后，检测河流交汇点处的拓扑关系并进行拓扑重构，从而保持河流交汇处的拓扑一致性。

D-P算法的原理是不断地构造新基线，计算点到基线的距离获得最大点线距的点并将其与预先设定的阈值比较，从而决定该点是否删除，以达到对曲线形态的化简^[6]。该方法的本质是通过点到基线的距离来表达点的形态特征。同样，面对三维曲线，构造的基线则应该是空间两点所构的空间直线段，而点到基线的距离则应是点到空间直线段的距离。如图 2所示，曲面ABCDEFGH为地形表面，曲线P₀P₁为河流曲线，直线P₀P₁为连接P₀、P₁的空间直线，平面ADEH为包含直线P₀P₁的平均坡面，P_i为曲线P₀P₁上任意一点，O_i为P_i垂直于直线P₀P₁的垂足，Q_i为P_i在平面ABCD上的投影点。显然，ΔP_iQ_iO_i为垂直于平面ADEH的直角三角形，设d_PO表示P_i、O_i两点的距离，则 ${{d}_{PO}}=\sqrt{{{({{d}_{PQ}})}^{2}}+{{({{d}_{QO}})}^{2}}}$ 。式中，d_PQ为P_i、Q_i两点的距离，表达P_i点在曲线P₀P₁中的高程形态；d_QO为Q_i、O_i两点的距离，表达P_i点在曲线P₀P₁中的二维形态；d_PO则综合了d_PQ和d_QO两者的含义，表达P_i点在曲线P₀P₁中的三维形态。因此，采用三维点线距扩展D-P算法能从曲线三维特征的角度选点，该算法的执行过程和D-P算法类似，此处不再赘述。

图 2 三维扩展DＧP算法的三维点线距表达 Fig. 2 Three-dimensional point-line distance principle of 3D extended D-P algorithm

图选项

水系树结构的构建关键在于建立干流与支流的关系，文献^[22]给出完整的河系树结构建立方法，但其构建的河系树结构只包含河流编号之间的父子关系，缺少河流交汇点信息，且该方法通过人为指定确定主干流，不利于全自动处理，尤其是当水系数据中包含多个河系，多条主干流时。为了便于水系中河流交汇点的描述与处理，本文设计了一个河流信息结构RS，其包括河流本身标识BID、该河流汇入的河流标识TID(“父”河流或干流河流)以及河流汇点坐标SP 3个要素。该结构的设计是基于每一条河流仅汇入一条干流河流或不汇入任何河流，如果是后者，则该河流为主干流。从树的逻辑结构来看，水系树结构是由节点组成的，每个节点包括一个数据域和若干个指针域，RS即为该节点的数据域，而其指针域则存储干流与支流的关系。

水系树结构的建立过程分为河流信息结构构造(RS构造)、主干流确定，干流-支流关系建立3步：①RS构造：RS包括BID、TID、SP 3要素，其中BID可直接从当前河流数据中获取，SP可通过比较河流两个端点高程的大小确定(高程值小即为河流汇点SP)，TID则通过遍历其他所有河流的所有点，找到与SP坐标一致的点，该点所在的河流标识即为TID；②主干流确定：由于主干流不汇入任何河流，因此当某河流RS结构中的TID数据为空(其值可设为-1)，则为主干流，故RS构造后主干流即可确定；③干流-支流关系建立：采用从上至下的层次构建方法，即从主干流开始，依次建立各干流-支流的关系。首先从主干流的RS结构中获取其河流标识，记为BID₀；然后依次遍历其他所有河流RS结构中的干流河流标识，记为TID_i；最后查找与BID₀相等的TID_i，将相匹配的河流记为主干流的一级支流。同理，根据一级支流确定二级支流，依次进行，整个水系树结构建立成功。图 3为水系示意图，其中数字1—10为各条河流，字母A—J则为各个交汇点，图 4为根据图 3构建的水系树结构图，其中字母A—J表示交汇点SP，第1个数字表示结点当前河流的编号BID，第2个数字表示汇入河流编号TID，其中，-1表示河流不存在。

图 3 水系示意图 Fig. 3 Schematic diagram of water system

图选项

图 4 水系树结构 Fig. 4 Tree structure of water system

图选项

干流化简后，原有交汇点处拓扑关系可能破坏，因此需要进行拓扑关系变化检测与拓扑重构，拓扑重构过程又分为干流交汇点重新生成和支流汇点移位两个过程。

交汇点处拓扑关系是否变化，是在干流化简后，判断该交汇点是否仍然存在于干流曲线上。由于单条河流曲线上高程值是单调变化的，从河流源点到汇点高程逐渐降低，因此，可以从干流源点开始依次遍历曲线上点，判断该点的高程Z_i与交汇点高程Z_sp的关系：如果Z_i=Z_sp，则遍历结束，交汇点存在，拓扑关系没有变化；如果Z_i＜Z_sp，则遍历结束，交汇点不存在，拓扑关系发生变化，此时，为了便于后续的拓扑重构，记录高程为Z_i的点在曲线中的位置，该位置将是新交汇点生成与插入的位置。 1.4.2干流交汇点重新生成 当干流中的交汇点删除后，需要生成新的交汇点，新的交汇点应取支流汇点到干流的旧交汇点相邻两点(未被删除)所构空间直线的垂足。如图 5所示，干流R₁的交汇点Q_i删除后，从支流R₂的汇点P_n作空间直线Q_i-1Q_i+1的垂足O为新的交汇点。考虑到空间直线Q_i-1Q_i+1相对于支流前进的方向有可能后退，即P_nO与P_nP_n-1所成夹角θ小于90°(图 5(b)，此时，支流流向不合常理，因此从P_n点逆流向上依次取支流上的点，重新计算相应的垂足O与角度θ，直到θ大于90°为止(图 5(b)中O′点为新交汇点)。

图 5 拓扑重构过程 Fig. 5 Process of topology reconstruction

图选项

新交汇点O生成之后，就需将支流的汇点移位到该点位置，因此将河流R₂的SP结构中的汇点SP重置为O点；同时，修改R₂的河流数据，根据在上述干流交汇点重新生成过程获得到支流汇点的前一点P_i(图 5(a)中为P_n，5(b)中为P_n-1)，将P_i后续点删除，并在P_i之后插入O点。

试验过程采用两份1∶5万水系数据，数据1是丘陵地貌区域的水系，其等高线间距为10 m，河流条数28条(图 6(a))；数据2是黄土地貌区域的水系，其等高线间距为20 m，河流条数21条(图 6(b))。试验采用D-P化简方法^[6]、分段D-P化简方法^[21]和本方法对数据进行9种尺度的化简(压缩率为0.1到0.9)。考虑到原始河流数据通常为二维数据，必须先对水系数据进行三维化处理，水系三维化方法采用文献^[23]提出的根据等高线进行反距离权重内插与分段线性内插相结合的方法，对河流源点与汇点采用反距离权重插值计算高程，对其他内部点先计算河流与等高线的交点，然后以相邻交点对河流分段，之后采用流水线线性内插方法计算高程。此外，为了能得到固定的压缩率，采用“智能阻尼震荡法”^[24]，通过一个可变步长智能调整阈值，直到化简结果满足预设的压缩率为止。算法步骤如下，给定初始阈值dis、阈值调整的初始步长step、一个目标压缩率rat_o和目标压缩率的容差tol： (1) 采用阈值dis进行化简并计算压缩率rat_c，并令diff =rat_c-rat_o； (2) 如果|diff|≤tol，结束；否则，分以下4种情况分别处理返回步骤(1)：①如果diff＜0，且上一次diff＜0，则dis=dis+step；②如果diff＜0，且上一次diff＞0，则step=step/2，dis=dis+step；③如果diff＞0，且上一次diff＞0，则dis=dis-step；④如果diff＞0，且上一次diff＜0，则step=step/2，dis=dis-step。

图 6 水系试验数据 Fig. 6 Experimental data of water system

图选项

线化简算法评价的研究很多，但主要集中在位置精度和形态结构特征两个方面^[25]。文献^[26]提出应从几何精度、结构形态、空间关系和适用范围4个方面评估曲线化简算法，其中前两者分别对应于位置精度与形态结构特征。因此，考虑到水系化简的实际情况，分别从几何精度、结构形态与空间关系3个方面对本方法进行评估。

几何精度评估的指标很多，如矢量位移、面积偏差、长度变化、曲折度等，各指标都存在一定的相关性。本文采用水系二维长度变化、三维长度变化(化简后水系长度减去原始水系长度)为评价指标，将本方法与D-P法和分段D-P法对比分析。

图 7为数据1的水系长度变化对比图，水系二维总长度为35 973.6 m，三维总长度为37 086.9 m。从二维长度变化来看，3种算法在各个压缩率上结果都相差不大(图 7(a))，相比较而言，D-P算法的变化最小；从三维长度变化来看，D-P算法和分段D-P算法相差不大，而本方法的变化则要远小于另外两种方法，且随着压缩率的增加差距更大(图 7(b))

图 8为数据2的水系长度变化对比图，河流二维总长度为87 109.1 m，三维总长度为87 433.6 m。从二维长度变化来看，3种算法差别很小，但仍然表现为随压缩率增大差距逐渐增大，相比较而言，本方法变化最小，D-P法次之(图 8(a))；从三维长度变化来看，本方法在各个压缩率上的变化要明显小于其他两种算法，D-P法的变化比分段D-P法略小(图 8(b))

图 7 3种化简方法对数据1的水系长度变化对比 Fig. 7 Contrast among length change of water system using three methods in data 1

图选项

图 8 3种化简方法对数据2的水系长度变化对比 Fig. 8 Contrast among length change of water system using three methods in data 2

图选项

评价曲线的形态结构特征可从分维数和特征点保持两方面进行^[26]。文献^[26]提出基于分形理论评估化简后曲线形态的一般方法，通过比较化简后曲线分维数与原始曲线分维数的相似程度，从而判定算法优劣。本文采用二维分维数与三维分维数为评价指标，前者基于二维曲线计算，后者基于三维曲线计算。分维数采用“脚规”法计算，通过不同步长量测曲线长度，并以线性回归模型拟合^{[17, 26]}。

表 1为本方法与D-P法化简后水系分维数计算结果，其中数据1的原始二维分维数为1.038 16，三维分维数为1.037 89；数据2的原始二维分维数为1.036 84，三维分维数为1.036 83。图 8为根据表 1生成的分维数变化(原始分维数减去化简后分维数)对比图。对于二维分维数，本方法与D-P法相差不大，仅在个别压缩尺度上，D-P法的变化更小；对于三维分维数，本方法在各个化简尺度数据上的变化都明显小于D-P法，尤其是数据1(图 9(a))表现更突出。

特征点保持主要考察不同化简方法对二维、三维特征点的保留情况。图 10为本方法与D-P法在同一化简尺度下水系特征点保留对比图(部分)，其中圆点为本方法保留特征点，三角形点为D-P法保留特征点。从二维特征点保持来看，两种方法差异不大，其保留的特征点都能较好呈现曲线的二维特征；然而，从三维特征点保持来看(借助地形)，本方法保留的特征点更好，尤其是像A、B等“叠水”位置，本方法保留的两点为曲线上三维拐点，显然优于D-P保留的两点。

图 9 DGP法与本方法化简后水系分维数对比 Fig. 9 Contrast between fractal dimension of water system using this method and D-P algorithm

图选项

图 10 本方法与DGP法化简后特征点保留对比图(部分) Fig. 10 Contrast between feature point of water system using this method and D-P algorithm

图选项

图 11 3种算法化简后拓扑一致性结果对比(部分) Fig. 11 Contrast among topological consistency of water system using three methods

图选项

本文提出了顾及拓扑一致性的三维水系化简方法，通过构建水系树结构、三维曲线扩展的D-P算法以及拓扑一致性检测与拓扑重构处理方法，实现了保持拓扑一致性的水系三维形态化简并进行了试验验证。结果显示：①该方法能保证河流间的拓扑关系，避免常规河流化简出现的拓扑不一致现象，且其拓扑一致性处理更合理；②该方法化简精度可靠，能最大限度地保证三维形态，随着综合尺度的增大，效果相对更优。该方法的不足之处在于常规河流数据大多不包含高程信息，需在化简之前先将河流曲线数据三维化，因此，设计一种快速高精度的水系高程信息提取方法非常关键，这也是后续的研究方向之一。