CN116431964A - 一种复杂河网水系骨架线生成的游程剥离法 - Google Patents

Abstract

本发明涉及一种复杂河网水系骨架线生成的游程剥离法；目的是克服现有技术内存占用大、计算效率低的缺陷。技术方案是：一种复杂河网水系骨架线生成的游程剥离法，步骤是：1、构建一个空的栅格场，其中每个栅格行生成一个空的游程链表；2、在每个行上完成栅格压缩，将属于河流范围的连续网格压缩成游程单元后插入栅格行游程链表；3、建立栅格场河流边界网格搜索模板；4、设置两个判别函数，判断栅格场上某网格是否是冗余网格；5、沿着栅格场边界搜索所有边界网格，分别基于两个判别函数标记并且剥离所有的冗余网格；6、分析是否有边界网格被剥离，若有，返回继续执行剥离，否则进入下一步；7、将未被剥离的网格转换为矢量形式的骨架线数据。

Description

一种复杂河网水系骨架线生成的游程剥离法

技术领域

本发明涉及一种水利信息技术，具体是一种复杂河网水系骨架线生成的游程剥离法。

背景技术

在水文学与水资源研究领域，河流系统是最基础也是最重要的研究对象。一般而言，水系通常是指流域内由干流和若干支流相互连接构成枝杈结构的河流，而河网则是指平原地区由多条河流相互连接构成网状结构的河流，水系和河网共同构成了河流系统。河流作为一种面状地物要素，其骨架线是对河流形状的抽象描述，通常指位于河流中心的曲线，它反映了河流的主延伸方向和形状特征。提取河流系统的骨架线有诸多应用场景，如在地图制图中经常要将双线河流化简为单线河流，在基于遥感影像的河流检测与识别领域也有河流骨架线提取的需求。同时，河流系统的骨架线数据可为水资源监测、土地规划、船舶运输等行业提供参考，具有重要的应用意义。

常见的骨架线提取方法主要分为矢量方法和栅格方法两类。其中矢量方法主要有Delaunay三角网法和Voronoi图法，栅格方法则有经典的细化算法以及数学形态学方法等。在矢量方法中，Delaunay三角网法能得到近似的骨架线，但是在处理类似河流等具有平行轮廓的多边形时，得到的骨架线会存在较多的“锯齿”；Voronoi图法不受噪声干扰，但算法设计较为复杂，且生成的骨架线需要“剪枝”。在栅格方法中，经典的栅格细化算法是一种特殊的多次迭代收缩算法，该算法在保持多边形拓扑不变性的条件下，逐步移除边界点或减少不必要点，直到推进至多边形内部提取出骨架线；数学形态学方法通过结构元素的组合对二值图像进行细化，具有算法简单、易于实现等优点，但是结构元素的选择需要根据栅格数据的特点设计，无法保证其普适性。总体而言，矢量方法计算效率较高，且更易于与面向目标的地图数据模型进行无缝结合，但是其程序设计较为复杂，很难避免一些后处理操作。与矢量方法相比，栅格方法能够较准确地获取多边形中轴，避免产生过多“毛刺”，但由于栅格方法需要事先对图形进行栅格化等操作，这在一定程度上降低了算法计算效率。

无论是栅格方法还是矢量方法，在应对诸如道路网络和河流系统时，如何提取那些分支众多、又长又窄的多边形的骨架线一直是研究的难点所在。特别是河流系统具有更为复杂的特性，以平原河网为例，大型河流的宽度可以达到几百米，而人工河道互相交织，最窄处宽度甚至只有不到5米，在应用栅格方法提取骨架线时必须慎重考虑栅格分辨率大小，需要尽可能加密网格以提升骨架线提取的精度。但是，更高分辨率的网格意味着算法对计算机内存使用和计算效率有着更高的要求。以一个100km×100km区域的河流系统为例，若网格分辨率设置为0.5米，按照每个网格存储一个字节的数据量计算，则栅格数据量将达到37.25GB，这对于个人计算机来说是一个严重的负担。一个更为典型的案例是对中国长江流域的水系提取骨架线，由于长江流域自东向西达到6000km，南北方向也达到了4000km，即使按照10米分辨率栅格来计算，其数据量也将达到223GB。

发明内容

本发明的目的是克服上述背景技术的不足，提供一种复杂河网水系骨架线生成的游程剥离法，该方法能够完全基于游程编码数据提取河流水系骨架线，具有内存占用少、计算效率高的特点。

本发明提供的技术方案是：

一种复杂河网水系骨架线生成的游程剥离法，包括如下步骤：

第一步、基于河流水系二值栅格数据，构建一个空的栅格场，设定空栅格场的参数；为每个栅格行生成一个空的游程链表。

第二步、按行读取河流水系二值栅格数据，并在行上完成栅格压缩，将属于河流范围的连续网格压缩成游程单元，并插入栅格行游程链表；

第三步、基于游程链表数据，考虑河流数据存在岛型结构复杂条件，建立栅格场河流边界网格搜索模板；

第四步、设置判别函数一、判别函数二，判断栅格场上某网格是否是冗余网格；冗余网格是指该网格被剥离后(网格值从1变成0)不会影响河流栅格数据的连通性；两个判别函数分别为判别函数一和判别函数二；

第五步、对游程链表数据，沿着栅格场边界搜索所有边界网格，基于判别函数一标记所有可以剥离的冗余网格，在搜索完毕后一次性剥离所有冗余网格；

第六步、对游程链表数据，沿着栅格场边界搜索所有边界网格，基于判别函数二标记所有可以剥离的冗余网格，在搜索完毕后一次性剥离所有冗余网格；

第七步、分析第五步和第六步是否有边界网格被剥离，若有，返回继续执行第五步，否则循环终止，进入第八步；

第八步，保留的未被剥离的网格即为河流骨架线栅格数据，通过栅格网格搜索，将栅格形式的骨架线数据转换为矢量形式的骨架线数据。

所述第一步中的河流水系二值栅格数据，是基于遥感数据提取得到，或者是通过河流矢量面数据栅格化得到。

所述空栅格场的参数为栅格范围、栅格行列数、网格宽和高参数，各参数值与河流水系二值栅格数据的栅格参数一致。

所述第一步中的游程链表中用于存储游程单元，每个游程单元都标记栅格行上游程单元跨越的起始行值和终止列值。

所述第二步中河流水系二值栅格数据，其中的“二值”指网格值为0或者1，其中0表示网格是背景网格，1表示该网格属于河流。

对于任意栅格行，栅格压缩和游程单元的生成步骤为：

(1)按照从左至右在栅格行上逐个网格遍历，查找网格的值，若是0则跳过，若发现是1，则将出现1的网格列值标记为起始列；

(2)继续按照从左至右搜索栅格行上的网格，一旦出现网格列值不为1，则将当前网格左侧值为1的网格列作为终止列；

(3)将起始列和终止列进行配对，生成一个游程单元，并插入游程链表，继续执行步骤(1)，直到栅格行上所有网格遍历完毕。

所述第三步中建立栅格场河流边界网格搜索模板，是采用一个2×2窗体对边界网格进行匹配，2×2窗体代表栅格场中四个网格，其中一个网格为当前搜索网格，其余三个网格为待搜索网格，根据四个网格的网格值分布，分析网格边界搜索的上、下、左、右走向，共存在24种河流边界网格搜索模板，其中内边界和外边界网格搜索模板分别为12个。

所述第四步中设置的两个判别函数，是采用一个3×3窗体对当前边界网格进行匹配，当前边界网格位于3×3窗体中间，基于判别函数判别当前网格是否可剥离，即是否可将网格值从1改成0而不影响河流栅格连通性。从当前边界网格正上方的网格开始，按照顺时针方向对相邻8个网格进行编号。对于8个网格中的任意一个网格，设置网格的0、1变换值，当某相邻网格值为0且其沿顺时针的下一网格值为1，则认为当前网格0、1变换值为1，否则为0。如此统计所有8个网格的0、1变换值。

所述判别函数一为：

式中，m和n分别为栅格场的总行数和总列数；

判别函数一所述的判别规则为：(1)、边界网格值为1，且8个相邻网格的0、1变换值之和为1；(2)、8个相邻网格值的和大于等于2，且小于等于6；(3)、1号网格、3号网格、5号网格的网格值之积为0；(4)、3号网格、5号网格、7号网格的网格值之积为0；

所述判别函数二为：

式中，m和n分别为栅格场的总行数和总列数；

判别函数二所述的判别规则为：(1)、边界网格值为1，且8个相邻网格的0、1变换值之和为1；(2)、8个相邻网格值的和大于等于2，且小于等于6；(3)、1号网格、3号网格、7号网格的网格值之积为0；(4)、1号网格、5号网格、7号网格的网格值之积为0；

采用D1和D2的判定每次交替进行。

所述第五步中沿着栅格场边界搜索所有边界网格，其具体步骤为：

(1)将栅格场中，所有游程单元的左右两侧进行标记，标记其是否有被搜索过，在初始状态下每个游程单元的左右两侧都处于未搜索状态；

(2)按照栅格行进行逐行遍历，查找栅格行中每个游程单元的左侧和右侧，若某游程单元左侧或右侧没有被搜索过，则将该游程单元左侧或右侧网格作为起始搜索网格，追踪方向为向上，同时将网格标记为已搜索状态；

(3)按照步骤三提供的边界搜索模板，执行网格边界搜索，查找搜索边界网格，并按照判别函数一标记所有边界网格是否需要进行剥离；

(4)完成边界网格搜索后，一次性将所有可以剥离的边界网格执行剥离操作。对于游程单元而言，其网格剥离与二值栅格数据不同的点在于，其不是将网格值从1变成0，是通过直接修改游程单元左侧或右侧的列值实现的；

所述第六步与第五步基本相同，区别仅在于判断网格是否可以剥离时，采用的是判别函数二。

所述第七步在于判断步骤五和步骤六是否能够继续剥离冗余的网格，如果无法剥离任何一个冗余网格，表示保留下的所有网格都位于河流骨架线上，此时必须终止步骤五和步骤六的循环判定操作。

所述第八步河流骨架线矢量数据的生成，步骤如下：

(1)将栅格场中保留的网格点分为三种：端点、连接点和分叉点；对所有网格点的类型进行标记；其中：

端点表示其8个相邻网格中，只有一个网格值为1；

连接点表示其8个相邻网格中，只有两个网格值为1，且这两个网格不相邻；或者有三个网格值为1，其中有两个网格处于相邻关系；

分叉点则表示其8个相邻网格中，有三个网格值为1，且这三个网格互不相邻；

(2)查找所有的属于端点的网格，记录网格中心点坐标，按照上、下、左、右、左上、左下、右下、右上的顺序搜索相邻网格，一旦找到未被搜索过的值为1的相邻网格，就将原网格标记为已搜索，记录新的网格中心点坐标，同步搜索相邻网格。循环执行上述步骤，直到搜索到另一个端点网格或分叉点网格，停止当前搜索；

(3)上述步骤完成后，所有端点处的网格已经完成顺序搜索并生成了骨架线，还有少量两端为分叉点的骨架线为完成搜索。随后查找所有未被搜索过的分叉点网格，按照步骤(2)的方式进行网格搜索并记录网格中心点坐标；

(4)将存储的所有骨架线坐标自动连接成线，并对子骨架线进行拼接，生成河流水系的整体骨架线。

本发明的有益效果是：

本发明针对当前栅格图像骨架线自动提取方法普遍存在的内存占用过高、计算效率低下的弊端，完全基于游程编码数据提取河流水系骨架线。该方法尝试将河流数据从直接栅格编码替换成游程编码存储，可以极大地压缩栅格数据量；采用游程边界追踪策略代替全栅格矩阵遍历，实现冗余像素的快速判断，从而能够大幅提升计算效率。本方法在内存耗用上不到经典栅格编码方法的1％，而计算效率是其10倍以上，完全具备在个人电脑上实现TB级河流栅格图像提取骨架线的应用潜力。

本发明的特点:

(1)本发明与常规直接栅格编码方法的优点是特别明显的，直接栅格编码下采用矩阵存储所有栅格值，特别是对于河流未覆盖的背景区域也存储了无效数值0，因此其数据量巨大。而本发明只用存在栅格行上河流覆盖的那些网格，且每个游程单元只存储位于河流边界处的网格值，可以极大的降低数据量。

(2)对于任一像素而言，由于无法知道该像素与周围八个像素的储值关系，因此现有的直接栅格编码方法，必须对整个二值图像进行循环遍历，按行逐个扫描网格；该方法效率较低，特别是对于那些值为0的背景像素，即使河流未覆盖该区域，仍无法避免判断该像素的值。而本发明因为游程编码数据只存储河流覆盖的网格，可以避免对那些背景像素进行遍历和开展冗余像素判别，进而可以极大地提升计算效率。

(3)现有的栅格细化方法是从外侧向内侧逐步剥离像素的，在直接栅格编码下，一次循环遍历时完全位于内部的像素是不会被剥离的。因此，直接栅格编码数据中存在大量的无效判断。而本发明采用游程编码数据，所有的网格遍历和判断都位于河流边界上，河流内部并不参与判断，这样可以在前述有益效果(2)的基础上进一步提升计算效率。

附图说明

图1是本发明实施例的流程示意图。

图2是某河道二值栅格数据示例；其中，b图是a图的局部放大图，c图是b图的局部放大图。

图3是执行游程剥离后栅格场中保留的河流骨架线网格示例；其中：a图是河流骨架线网格像素点，b图是河流骨架线网格像素点局部放大图，c图是河流骨架线网格像素点追踪后叠加了矢量形式骨架线的效果图。

图4是某地河流水系二值栅格数据骨架线提取效果图；其中：a图、c图、e图分别是局部放大的河流二值栅格图，b图、d图、f图分别是局部放大的河流二值栅格骨架线矢量数据图。

图5是现有的直接栅格编码数据及其冗余像素判定示意图。

图6是本发明的游程编码数据及其冗余像素判定示意图。

图7是某河流外边界和内边界示意图。

图8是本发明所述的游程边界搜索模板。

图9是本发明所述的游程边界搜索及其冗余像素剥离示意图。

图10是本发明所述河流骨架线栅格追踪端点、连接点和分叉点标记示意图。

图11是本发明所述的端点、连接点和分叉点追踪走向模板。

具体实施方式

为了能在普通个人计算机上实现大范围、高精度、高效率的河流系统骨架线自动提取，必须对经典栅格方法进行改进。考虑到河流数据在矢量转换栅格时是采用全栅格编码存储的，而这种二值栅格图存在大量的数据冗余。为此，本发明提出使用游程编码对河流数据进行存储，并基于游程编码数据实现河流骨架线的提取，可以将其称为游程剥离法。在直接栅格编码方式下，对于任意一行栅格，无论该行是否有河流跨越，都必须按照列号逐个存储像素值，而在游程编码方式下，只用考虑河流覆盖的那些网格，且仅仅需要标记每一栅格行中河流覆盖的起始列和终止列，这可以大大压缩数据存储量，进而实现河流骨架线的高精度、高效率生成。

下面结合附图，对本发明的技术方案进行清楚、完整地描述。

图1所示为本发明其中一个实施例的流程示意图。图2是某河道二值栅格数据示例，图中的任意一点代表的是栅格场中网格中心点的坐标，一般而言也可以将栅格场中的网格称为“像素”或“像素点”。图3为执行游程剥离后栅格场中保留的河流骨架线网格示例，可以看到执行游程剥离以后，栅格场上只保留了河流骨架线处的像素，其中c图中叠加了矢量形式的骨架线。图4为某地河流水系二值栅格数据骨架线提取效果图(其中：b图是a图的提取效果图，d图是c图的提取效果图、f图是e图的提取效果图)，可以从局部放大图中清晰地看到骨架线比原始栅格图要更细。

图1所示的本发明实施例复杂河网水系骨架线生成的游程剥离法，包括如下步骤：

第一步、基于河流水系二值栅格数据，构建一个空的栅格场，设定空栅格场的栅格范围、栅格行列数、网格宽和高参数；为每个栅格行生成一个空的游程链表。

此处，空的栅格场实际是用于游程编码数据表达河流栅格数据的一组参数，这些参数共同决定了栅格场的大小，具体参数包括栅格范围、栅格行列数、网格宽和高参数；其中栅格范围标记了栅格场的左下角(最小)和右上角(最大)坐标，栅格行列数决定了栅格场的网格密度，而网格宽和高参数由栅格场范围和栅格行列数决定，一般而言网格的宽和高相同。

如图2所示，该河流二值栅格数据中，河流覆盖的范围中网格值为1,未被河流覆盖的区域网格值为0，对于值为0的网格，其是作为背景网格存在。游程编码所代表的栅格场与当前河流二值栅格数据的栅格场参数完全相同。

游程单元和行游程链表的C语言形式如下所示，其中游程单元标记了改行游程跨越的起始网格列值和终止网格列值，游程链表中存储了游程单元数量，以及第一个游程和最后一个游程的内存地址。

第二步、按行读取河流水系二值栅格数据，并在行上完成栅格压缩，将属于河流范围的连续网格压缩成游程单元，并插入栅格行游程链表；

对于任意栅格行，栅格压缩和游程单元的生成步骤为：

(1)按照从左至右在栅格行上逐个网格遍历，查找网格的值，若是0则跳过，若发现是1，则将出现1的网格列值标记为起始列；

(2)继续按照从左至右搜索栅格行上的网格，一旦出现网格列值不为1，则将当前网格左侧值为1的网格列作为终止列；

(3)将起始列和终止列进行配对，生成一个游程单元，并插入游程链表，继续执行步骤(1)，直到栅格行上所有网格遍历完毕。

图5和图6分别是直接栅格编码以及游程编码存储示例(图中数字l1至l9代表网格行数)；在图5中每个栅格场上每个网格的值都进行了存储，现有的直接栅格编码是以矩阵的形式存储网格值的；而图6中游程编码数据是一种压缩形式，从图中可以看出每个栅格行上都只有一个游程单元(以长条矩形表达)，每个游程单元只用存储长条矩形覆盖的第一个网格的列值以及最后一个网格的列值即可。图5需要存储9×15＝135个值，而图6只用存储2×7＝14个值。可以看出，对于二值栅格数据，游程编码数据有极大的压缩比率，其不是存储的像素值，而是网格列值。

第三步、基于游程链表数据，充分考虑河流数据存在岛型结构等复杂条件，建立栅格场河流边界网格搜索模板。

图7是一个具有岛型结构的河流示意图，河流的内外边界共同组成了游程单元序列，要完整的从游程链表中搜索河流内外边界，必须建立边界网格搜索模板。建立栅格场河流边界网格搜索模板(如图8所示)，是采用一个2×2窗体对边界网格进行匹配，2×2窗体代表栅格场中四个网格，其中一个网格为当前搜索网格，其余三个网格为待搜索网格，根据四个网格的网格值分布，分析网格边界搜索的上、下、左、右走向，共存在24种河流边界网格搜索模板，其中内边界和外边界网格搜索模板分别为12个。

图9是基于游程数据的边界网格搜索示例；从图中可以看出，a图和b图分别是图8中第2行的“自上向上”和“自上向右”搜索，其中b图是a图搜索的下一步走向，即网格边界搜索方向从“自上向上”和“自上向右”搜索，此过程同步完成了边界网格的标记。

第四步、设置两个判别函数，用于判断栅格场上某网格是否是冗余网格，冗余网格是指该网格被剥离后(网格值从1变成0)不会影响河流栅格数据的连通性。两个判别函数分别称为判别函数一和判别函数二；

图5是现有的栅格细化方法示意图；在二值图像中将目标像素标记为1，背景像素标记为0。设p0是栅格场中位于x列y行的像素f(x，y)，以p0为中心像素，按照顺时针标记与p0相邻的八个像素点，分别为p1＝f(x，y+1)、p2＝f(x+1，y+1)、p3＝f(x+1，y)、p4＝f(x+1，y-1)、p5＝f(x，y-1)、p6＝f(x-1，y-1)、p7＝f(x-1，y)、p8＝f(x-1，y+1)。

定义任意两个相邻像素的基本关系式如下：

式中w为1至8中的某个值，比如设置w为1，则当p1＝0，且p2＝1时，

否则让/>

特别的，当w为8，因为没有p9，此时应该比较p8和p1的关系，即当p8＝0且p1＝1时，/>

统计像素p0的所有两两相邻像素关系式之和如下：

如公式(3)和(4)所示，J1(x，y)和J2(x，y)是针对栅格场上列为x行为y的像素点p0的两个判别函数，用来判断当前网格是否是冗余网格。在不满足限定条件时，J1(x，y)和J2(x，y)的值为0。若p0＝1，且满足公式(3)中J1(x，y)＝1或者公式(4)中J2(x，y)＝1，则p0是可删除点，即可以将p0设置为背景像素，此时不会影响栅格图像的连通性。

第五步、对游程链表数据，沿着栅格场边界搜索所有边界网格，判别函数一标记所有游程单元左右两侧可以剥离的边界网格，在搜索完毕后一次性剥离所有可剥离的边界网格；

按照公式(3)，遍历栅格场中所有像素，一次性删除满足J1(x，y)＝1的像素点。栅格场中像素连通性检测如公式(5)所示：

式中，m和n分别为栅格场的总行数和总列数，D1为遍历一次栅格场后，需要删除的像素点数量。根据上述描述，第一次边界搜索并删除冗余像素的伪代码如下所示：

If p1+p2+p3+p4+p5+p6+p7+p8≥2and p1+p2+p3+p4+p5+p6+p7+p8≤6；

Int Flag＝0；

If p1＝＝0and p2＝＝1Then Flag++；

If p2＝＝0and p3＝＝1Then Flag++；

If p3＝＝0and p4＝＝1Then Flag++；

If p4＝＝0and p5＝＝1Then Flag++；

If p5＝＝0and p6＝＝1Then Flag++；

If p6＝＝0and p7＝＝1Then Flag++；

If p7＝＝0and p8＝＝1Then Flag++；

If p8＝＝0and p1＝＝1Then Flag++；

If Flag＝＝1and p1*p3*p5＝＝0and p3*p5*p7＝＝0Then p0 can removed

第六步、对游程链表数据，沿着栅格场边界搜索所有边界网格，判别函数二标记所有游程单元左右两侧可以剥离的边界网格，在搜索完毕后一次性剥离所有可剥离的边界网格；

按照公式(4)，遍历栅格场中所有像素，一次性删除满足J2(x，y)＝1的像素点。栅格场中像素连通性检测如(6)所示：

式中，m和n分别为栅格场的总行数和总列数，D2为遍历一次栅格场后，需要删除的像素点数量。根据上述描述，结合公式(1)～(6)，栅格细化核心算法的伪代码如下所示：

If p1+p2+p3+p4+p5+p6+p7+p8≥2and p1+p2+p3+p4+p5+p6+p7+p8≤6；

Int Flag＝0；

If p1＝＝0and p2＝＝1Then Flag++；

If p2＝＝0and p3＝＝1Then Flag++；

If p3＝＝0and p4＝＝1Then Flag++；

If p4＝＝0and p5＝＝1Then Flag++；

If p5＝＝0and p6＝＝1Then Flag++；

If p6＝＝0and p7＝＝1Then Flag++；

If p7＝＝0and p8＝＝1Then Flag++；

If p8＝＝0and p1＝＝1Then Flag++；

If Flag＝＝1and p1*p3*p7＝＝0and p1*p5*p7＝＝0Then p0 can removed

第七步、分析第五步和第六步是否有边界网格被剥离，若有，返回继续执行第五步，否则循环终止，进入第八步；

在第五步和第六步中，采用D1和D2的判别每次交替进行；比如先执行D1，删除所有可以删除的像素点，再执行D2，同样删除所有可以删除的像素点，再次执行D1，如此循环，直到栅格场中已经没有可以删除的像素点。因此，第五步和第六步是配对出现的，一旦执行第五步，就一定会顺序执行第六步，不断往复循环。

第八步，保留的未被剥离的网格即为河流骨架线栅格数据，通过栅格网格搜索，将栅格形式的骨架线数据转换为矢量形式的骨架线数据。

如图10所示，在对游程单元两侧进行逐步剥离之后，最终保留下来的游程单元所代表的像素将出现如图中一样的效果。大量的游程单元只有一个像素点，即游程起始列值和终止列值相同。也有部分游程单元还保留有多个像素点，如分叉点所在行的游程。需要在这种游程数据中，设计栅格转换矢量策略，生成矢量形式的骨架线数据。

为此，本发明将这种骨架线所在的网格像素点分为三类，分别为端点、连接点和分叉点。其中，端点表示其8个相邻网格中，只有一个网格值为1，如图11中a图所示；连接点表示其8个相邻网格中，只有两个网格值为1，且这两个网格不相邻，或者有三个网格值为1，其中有两个网格处于相邻关系，如图11中b图所示；分叉点则表示其8个相邻网格中，有三个网格值为1，且这三个网格互不相邻，对所有的网格类型进行标记，如图11中c图所示。

图11是端点、连接点和分叉点追踪走向模板，其中九宫格中间的网格代表当前追踪网格像素点，8个相邻网格中网格中心为画点的为尚未追踪到网格，画点的表示该网格已经被追踪过，后续不允许再次追踪。

本发明的骨架线栅格像素点追踪策略为，首先从所有端点处开始追踪，如图11中a图所示，因为端点只有相邻网格像素点有值，如果该网格像素点未被追踪过则直接追踪到该网格像素点，若该网格已经被追踪过说明则直接退出此次追踪，将所有已经追踪的像素点顺序连接生成骨架线矢量点数据。同理，若追踪过程中发现追踪到如图11中c所示的分叉点，也直接退出此次追踪，将所有已经追踪的像素点顺序连接生成骨架线矢量点数据。

在从所有端点处完成追踪以后，栅格场中仍有少量骨架线未被追踪，这些骨架线的两个端点是分叉点。如图10所示，标记分叉点位置往下的一段骨架线即是此种情形。因此，从所有端点处完成追踪以后，需要再次搜索栅格场，找到那些未被追踪过的分叉点，将所有骨架线两端都是分叉点的网格单元追踪出来。最终，对追踪到的所有骨架线进行拼接，生成完整的河流骨架线矢量线数据。

Claims (10)

1.一种复杂河网水系骨架线生成的游程剥离法，包括如下步骤：

第一步、基于河流水系二值栅格数据，构建一个空的栅格场，设定空栅格场的参数；为每个栅格行生成一个空的游程链表；

第二步、按行读取河流水系二值栅格数据，并在行上完成栅格压缩，将属于河流范围的连续网格压缩成游程单元，并插入栅格行游程链表；

第三步、基于游程链表数据，考虑河流数据存在岛型结构复杂条件，建立栅格场河流边界网格搜索模板；

第四步、设置两个判别函数，判断栅格场上某网格是否是冗余网格；冗余网格是指该网格被剥离后不会影响河流栅格数据的连通性；两个判别函数分别为判别函数一和判别函数二；

第五步、对游程链表数据，沿着栅格场边界搜索所有边界网格，基于判别函数一标记所有可以剥离的冗余网格，并在搜索完毕后一次性剥离所有冗余网格；

第六步、对游程链表数据，沿着栅格场边界搜索所有边界网格，基于判别函数二标记所有可以剥离的冗余网格，在搜索完毕后一次性剥离所有冗余网格；

第七步、分析第五步和第六步是否有边界网格被剥离；若有，返回继续执行第五步，否则循环终止，进入第八步；

第八步，保留的未被剥离的网格即为河流骨架线栅格数据，通过栅格网格搜索，将栅格形式的骨架线数据转换为矢量形式的骨架线数据。

2.如权利要求1所述的复杂河网水系骨架线生成的游程剥离法，其特征在于：所述第一步中的河流水系二值栅格数据，是基于遥感数据提取得到，或者是通过河流矢量面数据栅格化得到；

所述空栅格场的参数为栅格范围、栅格行列数、网格宽和高参数，各参数值与河流水系二值栅格数据的栅格参数一致。

3.如权利要求2所述的复杂河网水系骨架线生成的游程剥离法，其特征在于：所述第一步中的游程链表用于存储游程单元，每个游程单元都标记栅格行上游程单元跨越的起始行值和终止列值。

4.如权利要求3所述的复杂河网水系骨架线生成的游程剥离法，其特征在于：所述第二步中的河流水系二值栅格数据，其中的“二值”指网格值为0或者1，其中0表示网格是背景网格，1表示该网格属于河流；

对于任意栅格行，栅格压缩和游程单元的生成步骤为：

(1)按照从左至右在栅格行上逐个网格遍历，查找网格的值，若是0则跳过，若发现是1，则将出现1的网格列值标记为起始列；

(2)继续按照从左至右搜索栅格行上的网格，一旦出现网格列值不为1，则将当前网格左侧值为1的网格列作为终止列；

(3)将起始列和终止列进行配对，生成一个游程单元，并插入游程链表，继续执行步骤(1)，直到栅格行上所有网格遍历完毕。

5.如权利要求4所述的复杂河网水系骨架线生成的游程剥离法，其特征在于：所述第三步中建立栅格场河流边界网格搜索模板，是采用一个2×2窗体对边界网格进行匹配，2×2窗体代表栅格场中四个网格，其中一个网格为当前搜索网格，其余三个网格为待搜索网格，根据四个网格的网格值分布，分析网格边界搜索的上、下、左、右走向，共存在24种河流边界网格搜索模板，其中内边界和外边界网格搜索模板分别为12个。

6.如权利要求5所述的复杂河网水系骨架线生成的游程剥离法，其特征在于：所述第四步中设置的两个判别函数，是采用一个3×3窗体对当前边界网格进行匹配，当前边界网格位于3×3窗体中间，基于判别函数判别当前网格是否可剥离，即是否可将网格值从1改成0而不影响河流栅格连通性；从当前边界网格正上方的网格开始，按照顺时针方向对相邻8个网格进行编号；对于8个网格中的任意一个网格，设置网格的0、1变换值；当某相邻网格值为0且其沿顺时针的下一网格值为1，则认为当前网格0、1变换值为1，否则为0；如此统计所有8个网格的0、1变换值。

7.如权利要求6所述的复杂河网水系骨架线生成的游程剥离法，其特征在于：所述判别函数一为：

式中，m和n分别为栅格场的总行数和总列数；

判别函数一所述的判别规则为：(1)、边界网格值为1，且8个相邻网格的0、1变换值之和为1；(2)、8个相邻网格值的和大于等于2，且小于等于6；(3)、1号网格、3号网格、5号网格的网格值之积为0；(4)、3号网格、5号网格、7号网格的网格值之积为0；

所述判别函数二为：

式中，m和n分别为栅格场的总行数和总列数；

判别函数二所述的判别规则为：(1)、边界网格值为1，且8个相邻网格的0、1变换值之和为1；(2)、8个相邻网格值的和大于等于2，且小于等于6；(3)、1号网格、3号网格、7号网格的网格值之积为0；(4)、1号网格、5号网格、7号网格的网格值之积为0；

采用D1和D2的判定每次交替进行。

8.如权利要求7所述的复杂河网水系骨架线生成的游程剥离法，其特征在于：所述第五步中沿着栅格场边界搜索所有边界网格，其具体步骤为：

(1)将栅格场中，所有游程单元的左右两侧进行标记，标记其是否有被搜索过，在初始状态下每个游程单元的左右两侧都处于未搜索状态；

(2)按照栅格行进行逐行遍历，查找栅格行中每个游程单元的左侧和右侧，若某游程单元左侧或右侧没有被搜索过，则将该游程单元左侧或右侧网格作为起始搜索网格，追踪方向为向上，同时将网格标记为已搜索状态；

(3)按照步骤三提供的边界搜索模板，执行网格边界搜索，查找搜索边界网格，并按照判别函数一标记所有边界网格是否需要进行剥离；

(4)完成边界网格搜索后，一次性将所有可以剥离的边界网格执行剥离操作；对于游程单元而言，其网格剥离是通过直接修改游程单元左侧或右侧的列值实现的；

所述第六步与第五步基本相同，区别仅在于判断网格是否可以剥离时，采用的是判别函数二。

9.如权利要求8所述的复杂河网水系骨架线生成的游程剥离法，其特征在于：所述第七步在于判断步骤五和步骤六是否能够继续剥离冗余的网格；如果无法剥离任何一个冗余网格，表示保留下的所有网格都位于河流骨架线上，此时必须终止步骤五和步骤六的循环判定操作。

10.如权利要求9所述的复杂河网水系骨架线生成的游程剥离法，其特征在于：所述第八步河流骨架线矢量数据的生成，步骤如下：

(1)将栅格场中保留的网格点分为端点、连接点和分叉点三种；对所有网格点的类型进行标记；其中：

端点表示其8个相邻网格中，只有一个网格值为1；

连接点表示其8个相邻网格中，只有两个网格值为1，且这两个网格不相邻；或者有三个网格值为1，其中有两个网格处于相邻关系；

分叉点则表示其8个相邻网格中，有三个网格值为1，且这三个网格互不相邻；

(2)查找所有的属于端点的网格，记录网格中心点坐标，按照上、下、左、右、左上、左下、右下、右上的顺序搜索相邻网格，一旦找到未被搜索过的值为1的相邻网格，就将原网格标记为已搜索，记录新的网格中心点坐标，同步搜索相邻网格；循环执行上述步骤，直到搜索到另一个端点网格或分叉点网格，停止当前搜索；

(3)上述步骤完成后，所有端点处的网格已经完成顺序搜索并生成了骨架线，还有少量两端为分叉点的骨架线为完成搜索；随后查找所有未被搜索过的分叉点网格，按照步骤(2)的方式进行网格搜索并记录网格中心点坐标；

(4)将存储的所有骨架线坐标自动连接成线，并对子骨架线进行拼接，生成河流水系的整体骨架线。

Images