一种基于TIN数据的河网提取方法

摘要

本发明涉及一种基于TIN数据的河网提取方法，属于地理信息与水文学技术领域。本法利用TIN三角形的梯度判断水流方向，在求得流向的基础上，通过判断TIN三角形每条边与水流流向向量的关系将空间三角形分为谷线三角形和非谷线三角形，并分别依据谷线三角形和非谷线三角形提取河网，采用从流域出口回溯的方法计算平坦区域三角形的流向以进一步提取出河道。本发明解决了通过判断TIN三角形每条边的汇水性来识别河道所定义的河道单一的问题以及无法提取平坦区域河道的问题。

权利要求

1.一种基于TIN数据的河网提取方法，其特征在于：包括如下步骤：步骤A：读取TIN数据，

步骤B：计算水流流向向量：

步骤B-1，利用梯度原理计算TIN中空间三角形重心处的梯度向量；

步骤B-2，定义以三角形重心为起点，沿与梯度向量方向相反的射线为水流流向向量，步骤C：提取谷线，将空间三角形分为谷线三角形和非谷线三角形，包括如下步骤：C-1，取TIN中任意一条边，在所取TIN边的左右三角形中，分别做所选边的方向向量与三角形水流流向向量的向量积；

当向量积小于零时，进一步计算所选边与水流流向向量的交点并根据所选边上点的取值范围判断交点是否在所选边的延长线上；

C-2，若左右两个三角形中所作向量积均小于零，所选边为河网谷线，含有河网谷线的三角形为谷线三角形；

若三角形的所有边均不是河网谷线，该三角形为非谷线三角形；

步骤D，提取河道，当步骤C的判别结果为非谷线三角形的，进入步骤D-1；当步骤C的判别结果为谷线三角形的，进入步骤D-2；

D-1，对于非谷线三角形，连接该非谷线三角形重心到其水流方向所指三角形的重心，构成重心—重心河道；

D-2，对于谷线三角形河道的提取，包括如下步骤：D-2-1：提取TIN中的所有谷线，并将谷线标记为谷线河道；

D-2-2：判断三角形包含的谷线数目，若有两条谷线，进入步骤D-2-2-1；若只有一条谷线，进入步骤D-2-2-2；

D-2-2-1：三角形有两条谷线，则这两条谷线与水流流向向量有且只有一个交点，连接三角形重心到该交点，构成重心—谷线河道，所述两条谷线不含延长线；

D-2-2-2：当水流流向向量与谷线相交时，连接三角形重心和交点，构成重心—谷线河道；当水流流向向量与谷线的延长线相交时，连接谷线上离交点较近的端点和三角形的重心，构成重心—谷线河道；

步骤E：提取平坦区域的河道，采用从流域出口回溯的方法计算平坦区域三角形的水流流向以进一步提取出河道。

说明书

一种基于TIN数据的河网提取方法

技术领域

[0001] 本发明涉及一种基于TIN数据的河网提取方法，属于地理信息与水文学技术领域。

背景技术

[0002] 随着空间技术和计算机技术的迅猛发展，水文研究也逐步迈向数字化，特别是3S技术(RS、GIS和GPS)的飞速发展为数字流域水系的研究提供了强有力技术支撑。遥感(RS)可以提供流域范围内不同时间、不同空间分辨率的格网数据，全球定位系统(GPS)能提供流域内水文实体(如河流、湖泊、水库等)精确空间位置信息，地理信息系统(GIS)能够采集、编辑和管理流域空间数据，此外还能进行空间分析、数据查询、地形分析等。RS、GIS和GPS等技术的相融合，可以为流域研究和管理提供自然地理、水文水资源、生态环境和社会经济等方面的快速、准确、实时、大范围观测的空间数据。3S技术相结合所提供的大量空间数据，使得水文研究，特别是数字流域的研究得到了进一步深入。随着水文研究的数字化，传统的河网已经不能满足水文研究的需要，数字河网概念应运而生。数字流域河网的提取已成为数字流域、分布式水文模型等水文研究的关键技术环节。

[0003] 目前关于数字河网的提取方法主要是针对格网DEM的，其主要有两类方法：第一类是：基于谷点的河网构建方法，先提取出格网中的谷点，再把这些谷点连接起来，从而得出河道和河网；该类方法简单，但是提取的河网存在不连续等问题，需要重新连接、修整才能生成合理的河网，故该方法运用较少。第二类是基于流向的构建方法，先确定每个格网的水流方向，再根据这些水流方向确定河网。基于谷点的河网提取。基于流向的河网提取方法的思路是，首先对研究区域进行填洼和平坦区域处理，再计算出每个DEM格网的水流流向，然后根据格网流向算出流向每个格网的汇水面积(上游汇水面的格网数目)，设定一个汇水面积阈值，将汇水面积大于或等于阈值的格网判定为流域河网的组成部分，最后按照水流方向连接这些格网形成最终的流域河网。由于格网DEM的结构简单，处理方便，基于格网DEM的河网提取算法研究较多，但是由于格网DEM的固有特性，利用格网DEM数据构建数字河网仍然存在较大的缺陷，特别是进行了洼地和平坦区域处理后，往往丢失了原来的地形信息，从而导致所构建的数字河网失真。相对于格网DEM而言，TIN结构的DEM数据具有良好的矢量性，能够更加精确地模拟地形地貌，因此，在河网提取的实际应用当中，TIN结构的DEM数据更具优势，可以利用三角形的梯度唯一地确定水流方向，河网可以作为三角形的特征模式来研究。Frank等于1986年首先提出基于TIN的流域河网提取，认为河网是由所有的在TIN上的谷线所组成的，除了谷线外，TIN还包含脊线(在TIN中认为是分水线)，在TIN中表现为三角形的边，其中谷线和脊线统称为特征线。卢振千等于2001年提出了一种融合等高线和规则格网有点的DEM-TIN模型，该方法要求规则格网内的河流必需位于由该格网生成的TIN的边上，并基于该模型进行了地形因子的提取(卢振千，吕妙儿，黄杏元.不规则三角网TIN在流域坡面汇流分析中的应用[J]测绘科学，2001，26(4)：45～48.)。Jones，N.L等于1990年研究了基于梯度的TIN的流向与边的汇水性计算算法，但是未进一步提取流域河网(Jones，N.L.，Wright，S.G.，and Maidment，D.R.Watershed delineation with triangle-based terrain models[J]，Journal of Hydraulic Engineering，1990，116：1232～1251.)。SidiYu等于1996年研究了边的汇水性，以及基于TIN的河网查询的算法复杂度，在SidiYu和Jones，N.L基础上，刘学军等于2008年提出了基于不规则三角网的河网水系提取算法，利用三角形的几何性质，提取出TIN的特征线，连接谷线形成河网。目前基于TIN的算法大多采用Frank的思想，认为河流在TIN的边上。

[0004] 上述算法均利用TIN三角形的梯度判断水流方向，在求得流向的基础上，通过判断TIN三角形每条边的汇水性来识别河道，河道的定义单一。当出现大面积平坦区域时，TIN数据中的三角形表现为导数为零，无法计算流向，进而不能提取河道。

发明内容

[0005] 本发明所要解决的技术问题是针对上述背景技术的不足，提供一种基于TIN数据的河网提取方法。采用从流域出口回溯的方法计算平坦区域三角形的流向以进一步提取出河道。

[0006] 本发明为实现上述发明目的采用如下技术方案：

[0007] 一种基于TIN数据的河网提取方法，包括如下步骤：

[0008] 步骤A：读取TIN数据，

[0009] 步骤B：计算水流流向向量，

[0010] 步骤C：提取谷线，将空间三角形分为谷线三角形和非谷线三角形，包括如下步骤：

[0011] C-1，取TIN中任意一条边，在所取TIN边的左右三角形中，分别做所选边的方向向量与三角形水流流向向量的向量积；

[0012] 当向量积小于零时，进一步计算所选边与水流流向向量的交点并根据所选边上点的取值范围判断交点是否在所选边的延长线上；

[0013] C-2，若左右两个三角形中所作向量积均小于零，所选边为河网谷线，含有河网谷线的三角形为谷线三角形；

[0014] 若三角形的所有边均不是河网谷线，该三角形为非谷线三角形；

[0015] 步骤D，提取河道，当步骤C的判别结果为非谷线三角形的，进入步骤D-1；当步骤C的判别结果为谷线三角形的，进入步骤D-2；

[0016] D-1，对于非谷线三角形，连接该非谷线三角形重心到其水流方向所指三角形的重心，构成重心-重心河道；

[0017] D-2，对于谷线三角形河道的提取，包括如下步骤：

[0018] D-2-1：提取TIN中的所有谷线，并将谷线标记为谷线河道；

[0019] D-2-2：判断三角形包含的谷线数目，若有两条谷线，进入步骤D-2-2-1；若只有一条谷线，进入步骤D-2-2-2；

[0020] D-2-2-1：三角形有两条谷线，则这两条谷线与水流流向向量有且只有一个交点，连接三角形重心到该交点，构成重心-谷线河道，所述两条谷线不含延长线；

[0021] D-2-2-2：当水流流向向量与谷线相交时，连接三角形重心和交点，构成重心-谷线河道；当水流流向向量与谷线的延长线相交时，连接谷线上离交点较近的端点和三角形的重心，构成重心-谷线河道；

[0022] 步骤E：提取平坦区域的河道，采用从流域出口回溯的方法计算平坦区域三角形的水流流向以进一步提取出河道。

[0023] 所述基于TIN数据的河网提取方法步骤B的具体实现如下：

[0024] 步骤B-1，利用梯度原理计算TIN中空间三角形重心处的梯度向量；

[0025] 步骤B-2，定义以三角形重心为起点，沿与梯度向量方向相反的的射线为水流流向向量。

[0026] 本发明采用上述技术方案，具有以下有益效果：

[0027] (1)解决了通过判断TIN三角形每条边的汇水性来识别河道所定义的河道单一的问题；

[0028] (2)解决了平坦区域河道提取的问题。

附图说明

[0029] 图1为基于TIN数据河网提取方法的流程图。

[0030] 图2为TIN的基本特征图。

[0031] 图3为标量场的梯度示意图。

[0032] 图4为三角形的梯度示意图。

[0033] 图5为空间三角形的流向向量示意图。

[0034] 图6为空间三角形边与流向的关系示意图。

[0035] 图7为空间三角流向向量与边的相交关系示意图。

[0036] 图8为相邻空间三角形公共边的分类示意图。

[0037] 图9为空间三角形边与相邻空间三角形的关系示意图。

[0038] 图10为谷线三角形的河道提取示意图。

[0039] 图11为实验测试数据示意图。

[0040] 图12为实验数据中研究区域的谷线分析结果图。

[0041] 图13为重心-重心河道实验结果图。

[0042] 图14为重心-谷线河道实验结果图。

[0043] 图15为回溯后的河道示意图。

具体实施方式

[0044] 下面结合附图对发明的技术方案进行详细说明：

[0045] 如图1所示，基于TIN数据河网提取方法的流程包括如下步骤：

[0046] (1)读取TIN数据；

[0047] (2)计算水流流向向量；

[0048] (3)提取谷线，将空间三角形分为谷线三角形和非谷线三角形；

[0049] (4)分别对谷线三角形和非谷线三角形进行河网提取；

[0050] (5)处理平坦区域，提取平坦区域的河道。

[0051] 如图2所示：在TIN模型中，有三种基本元素，分别为：点、边和面。点是数据采集时的数据采集点，在TIN结构中表现为相邻三角形的公共顶点，也是在进行河网提取时需要的数据基础。在TIN模型中点是具有三维空间信息的点，表现为(x，y，z)，其中z表示高程。边是TIN模型中的相邻三角形的公共边界，同时也是整个研究区域的边界，它把孤立的数据点关联起来，包含了大量地形特征。在数学上边是空间线段，不具备可导性质。根据TIN的就近原则，面是地形空间中三个相距最近的点构成三角形面，是用于数字化表达地形的基本单元，在进行河网提取时将用到其一重要的性质-梯度。

[0052] 梯度是一个重要的几何性质，在向量微积分学中，梯度是标量场的一个向量场，标量场中某点的梯度表现为经过该点的一个向量，沿着该向量的方向，标量增长的速度最快，如图3所示。对于二维平面区域上的函数z＝f(x，y)的梯度的数学定义为：

[0053] 设函数z＝f(x，y)在点(x，y)处的方向导数存在，则z＝f(x，y)在(x，y)处的梯度记为：

[0054]

[0055] 图4给出了TIN中的空间三角形中的点(x，y)处的梯度。它是以点(x，y)为起点的射线，其物理意义是沿此方向高程值z上升最快。在实际的地形中，如果高程值z用函数z＝f(x，y)表示，那么点(x，y)处的梯度就是坡度最陡峭的方向。在等高线模型中，由于等高线越密集，坡度越陡峭，故梯度总是指向等高线密集的。在TIN模型中，一个特定的空间三角形的梯度是一组平行的向量，它们可以由空间三角形的方程求出，这组向量的方向就是三角形所表达地形最陡峭的方向。

[0056] 本发明的方法是通过如下步骤实现的：

[0057] A.打开TIN数据，读入算法所需要的地形数据。

[0058] B.利用梯度原理计算TIN数据中空间三角形重心处的梯度，并进一步计算水流方向。具体包含如下步骤：

[0059] B-1利用TIN中空间三角形的顶点坐标计算出空间三角形的梯度：

[0060] 设空间三角形的三个顶点分别为：A(xa，ya，za)，B(xb，yb，zb)，C(xc，yc，zc)，则该三角形的梯度可以按如下步骤计算。

[0061]

[0062] 由于梯度的计算需要求偏导数，且数字高程模型的数学定义z＝f(x，y)，故对式2进行进一步化解得：

[0063] ax+by+cz+d＝0(3)

[0064] 其中

[0065] 当c≠0时，对式3进行变形，表达为z关于x，y的函数：

[0066]

[0067] 这样，就把XOY平面上的点坐标(x，y)映射到了空间高程z上面，在c≠0的条件下，z＝f(x，y)在其定义域内的方向导数存在，对式5分别求关于x，y偏导数可得：

[0068]

[0069]

[0070] 由导数的定义可以得到函数z＝f(x，y)的梯度：

[0071]

[0072] 令：

[0073] 式9中的向量 就是空间三角形平面内高程上升最快的方向，即最陡峭的方向。

[0074] B-2TIN中空间三角形流向的确定

[0075] 空间三角形的梯度向量，它表示该三角形表达地形的高程值沿着该方向上升最快，反之，式9向量方向的反方向便是高程值下降最快的方向。

[0076] 令： (c≠0) (10)

[0077] 式10给出了高程值下降最快方向的表达式，结合物理意义，定义式10中 为空间三角形的水流方向向量。由向量的定义可知，在三维空间内，式10其实是一组平行向量，称其为平行流向向量，如图5中虚线箭头所示。利用几何知识可知，在已知方向的情况下，要唯一确定一个向量，还需要该向量上的一点。由于本发明针对的是空间三角形的水流方向，水是在重力的作用下流动，故选择经过空间三角形ΔABC的重心点O(xo，yo，zo)的向量作为ΔABC的水流方向向量，称其为ΔABC的流向向量，如图5中实线箭头所示。特别地，为了更加精确地表示水流方向，把流向向量定义为以ΔABC重心为起点、方向为向量 的射线 记为水流方向。

[0078] C.谷线的提取

[0079] C-1水流方向与TIN空间三角形边向量关系的计算如下：

[0080] 为了减少算法的计算量，按照数学原理将TIN的空间三角形投影到二维平面，如图6所示。这样就可利用水流流向向量和边的方向向量来确定水流方向与边的关系。如附图6所示，设三角形顶点按逆时针方向排列，任意一边如AB边，起点和终点分别为A(xa，ya)和B(xb，yb)，其方向向量记为 方向向量为(xb-xa，yb-ya)，则AB边与该三角形水流方向的关系为：

[0081]

[0082] 如果IAB＞0，则水流离开该边；如果IAB＜0，则水流流向该边(水流向量可能与三角形边的延长线相交)；如果IAB＝0，则水流方向与该边方向相同(水流方向与该边方向向量平行)。

[0083] 由于在利用式11计算边与水流方向的关系时，把边当成是直线(利用的是直线的方向向量)，而不是线段，因此一般情况下，按式11计算，会出现两条边的关系值I小于0的情况：如图7所示，水流方向向量 与边AB和CA的I值均小于0，其中一个是水流方向向量交于边的延长线，与CA的延长线AA′相交。对于式12中IAB＜0情况下，为了能够区分出水流方向向量与边延长线相交的情况，需要把边的表示方式由直线改为线段，所采用的方法是：先求出水流方向向量与边(包括延长线)的交点，然后判断交点是否在边上(不包含延长线)。要确定地计算出水流方向向量与三角形的哪条边相交，还需要在式11的基础上作进一步地计算。

[0084] 如图7所示：设边AB的起点为A(xa，ya)，终点为B(xb，yb)，则边AB的方程为：

[0085] α1x+β1y+χ1＝0(12)

[0086] 其中：α1＝yb-ya，β1＝xa-xb，χ1＝ya×(xb-xa)-xa×(yb-ya)

[0087] 设水流方向向量 的起点为O(xo，yo)，方向向量 则其方程为：

[0088] αox+βoy+χo＝0(13)

[0089] 其中：αo＝b，βo＝a，χo＝a×yo-b×xo

[0090] 在IAB＜0的情况下，联立式12与式13，可得 与AB或其延长线的交点：

[0091]

[0092] xmin≤xi≤xmax (15)

[0093] ymin≤yi≤ymax (16)

[0094] 若式15和式16均成立，则 与AB相交，否则 与AB的延长线相交。

[0095] 在进行正式提取前需要对TIN数据中的空间三角形进行谷线分析，以为后续工作做准备。基于Frank的思想，本方法认为TIN中的谷线是河网的一部分，如图8所示，边AB是否为谷线，可以通过计算其左右三角形TL和TR的流向与边AB的关系来确定。

[0096] C-2设TIN中任意一条边为AB，其左右三角形分别为TL和TR，相应的水流方向和AB边的关系为IL和IR。

[0097] IL≤0，IR＜0或IL＜0，IR≤0 (17)

[0098] 若式17成立，则边AB为河网谷线，TL和TR为谷线三角形，如图8(a)所示；若IL＞0且IR＞0，则边AB则为河网脊线，如图8(b)所示；否则边AB位非特征线，如图8(c)所示。若三角形的所有边均不是河网谷线，即该三角形只含有河网脊线和非特征线，则该三角形为非谷线三角形。

[0099] 通过对谷线的提取可以将TIN中的三角形分为谷线三角形和非谷线三角形两类，这两类三角形的水流流向计算方法有所区别，形成的河道也不一样；前者的河道从其重心出发到其流出三角形的重心，后者的河道从三角形的重心出发沿水流方向到谷线结束。

[0100] D.TIN中三角形河道的提取

[0101] D-1非谷线三角形中河道的提取思想是：若三角形ΔABC的流向向量所指三角形为ΔABC′，这样连接ΔABC的重心O(xo，yo，zo)和ΔABC′的重心O′(x′o，y′o，z′o)，就形成了ΔABC和ΔABC′间的简单河道，记为重心-重心河道，如图9所示的OO′。

[0102] D-2对于谷线三角形，若直接连接重心形成河道，水将在这两个三角形间循环，针对这种情况，则选择重心到谷线的交点作为河道，记为重心-谷线河道，如图10(c)所示的OLIL和ORIR)。

[0103] 谷线三角形中河道提取的具体步骤如下：

[0104] D-2-1：提取TIN中的所有谷线。

[0105] D-2-2：判断三角形T包含的谷线数目。

[0106] D-2-2-1：若三角形T含有两条谷线Ei，Ej(如图10(a)所示)，则 与Ei，Ej有且只有一个交点PI，连接T的重心OT与PI形成河道；

[0107] D-2-2-2：若三角形T只含有一条谷线Ei(如图10(b)所示)，则进一步讨论 与Ei(含延长线)的相交关系；

[0108] ′若 与Ei相交(不含延长线，交点记为PI)连接T的重心OT与PI形成河道，若与Ei的延长线相交，设距离交点PI较近Ei的端点为PE，连接T的重心OT与PE形成河道。

[0109] 在三角形T只含有一条谷线Ei且 与Ei的延长线相交的情况下，由于 的方向是高程较低处，PE一定是Ei高程较低的端点，故连接OT与PE形成河道是合理的。

[0110] E平坦区域河道的提取

[0111] 由于平坦区域空间三角形的顶点高程值相等，三角形平面内所有点的导数均为零，从而无法进行水流方向的计算。本步骤在步骤C.2的基础上，采用从流域出口回溯的方法计算平坦区域三角形的流向以进一步提取出河道。具体思想是：从流域出口进行回溯，遍历流域出口点的相邻的三角形，若有方向没有确定的，则设置其流出三角形为流域出口三角形，依次回溯上去。比如取某三角形T，则遍历其相邻三角形，若其相邻三角形有三角形Tn的流出三角形没有确定，则将Tn的流出三角形设置为T，像这样直至所有三角形的方向都确定为止，以此确定平坦区域的河网河道。

[0112] 下面举一个具体例：我们选择了一组TIN地形数据进行试验，该组数据具有3659个数据点、6211个空间三角形共9905条边。分别设置汇水面积阈值500和2000个单位。

[0113] 本方案可以结合计算机系统实现。

[0114] 为了验证本方案的实际效果，我们做了如下试验：

[0115] 第一步：选取了一组如图11所示的试验数据，读入计算机采用算法定义的数据结构进行图11(a)为实验数据的地形图，图11(b)为实验数据的TIN结构。

[0116] 第二步：进行流向计算。

[0117] 第三步：谷线分析。谷线分析结果如图12所示(限于地形数据分辨率，仅显示图12中研究区的谷线分析结果)，图中白色线即为谷线。本方法能够提取TIN数据中的谷线，对于输入的TIN数据能够准确地提取出谷线，脊线和非特征线，且提取的谷线符合实际地形。

[0118] 第四步：河道提取。河道提取结果如图13所示，该结果仅提取重心-重心河道，造成很多重心-重心河道是分离且不连续的(河道AB)，图14是对谷线三角形进行河道提取后的效果图，可以发现，本方法在处理了谷线三角形后，形成的重心-谷线河道BE将重心-重心河道AB与谷线连接起来，从而将分离的河道相互连通。

[0119] 第五步：平坦区域的河道提取。如图15所示在对选定的实验数据进行回溯的过程中，处理了图15中平坦区域三角形T，在回溯到三角形T1时确定三角形T的流向，提取出重心-重心河道FG。分析图15可以发现本方法可以提取平坦区域的河道。

[0120] 从以上实验分析可以得出，本方法能够实现对TIN数据的流域河网提取，和现有的基于TIN数据的河网提取方法相比，本发明提出的方法能够实现基于TIN数据的河道提取，河道定义丰富，能够平坦区域的河道提取问题，而且在提取精确度上具有较大优势。