CN113989653B - 一种河流平面几何指标与拓扑结构的提取方法 - Google Patents

Abstract

本发明公开了一种河流平面几何指标与拓扑结构的提取方法，其IPC分类号为E02B1/00。该方法基于复杂网络分析技术，以卫星图片为输入量，由河道水体的提取、河流中心线骨架的刻画，河流分汊点与端点的识别与检测和河流网络拓扑结构与平面几何指标的提取四个主要部分组成，在考虑了江心滩，沙坝等对于河道平面形态所产生的影响的同时，尽可能保证了所提取的河道的连续性，因此对具有复杂平面形态的辫状河流，畸形河湾等同样具有较强的适用性。本发明不仅节约了现场测量所带来的大量人力物力成本，还兼顾了所获取数据的准确性与多样性，为河流研究提供了技术支持。

Description

一种河流平面几何指标与拓扑结构的提取方法

技术领域

本发明属于水利工程技术领域，具体涉及一种基于复杂网络分析技术的河流平面几何指标与拓扑结构的自动化提取方法，其IPC分类号为E02B1/00。

背景技术

地球上的河流系统通常具有不同的平面形态，这往往被认为是与河流行为密切相关的。河床演变学的主要内容之一就是研究上述两者之间的关系，为此需要定义包括河道宽度、河道长度、分汊角、弯曲度和迁移速率等在内的一系列量化指标来区分和描述不同的河流平面形态。因此，高效准确地绘制和评估河流平面形态并且从中提取出相应的几何指标，一直以来是河流研究的热点问题之一。

目前，传统的河流平面几何指标提取方法仍然依赖于野外测量，然而，这种消耗大量人力物力的方式无法保证河流数据在时空分布上的连续性；随着遥感技术(RemoteSensing,RS)与图像处理技术(Image Processing,IP)的发展，从河流高分辨率图像中获取平面几何指标的方法已经得到了广泛地应用，大量的自动化提取算法应运而生，例如，RivWidth、Rivamap、PyRIS、RiMARS等，然而，上述河流平面指标提取方法仍存在诸多不足：首先，他们很难识别河流中的江心滩与沙坝；其次，对于所检测出的河道，其连续性无法得到保证。这意味着上述河流平面几何指标自动化提取方法对于具有复杂平面形态的河型，例如辫状河流，并不适用。

与此同时，由于河流可以被看做是一种由节点(node)与链条(link)组成的网络结构，因此基于图论的复杂网络分析方法已经被引入到了河流平面指标的提取研究中。最近的研究已经实现了河流卫星图片中河网拓扑结构的构建，它虽然可以解决所提取河道的连续性问题，但目前仍停留在河网拓扑结构的提取与分析中，所提取的有效指标也只包含有河流中河道与分汊点的数目，尚未实现对于河流平面几何指标的自动化提取。因此，对于基于复杂网络分析的河流平面几何指标自动化方法的研究与提出，不仅可以节约大量的人力物力资源，同时也保证了所得数据的准确性与全面性，为河流研究提供了技术支撑。

发明内容

本发明的目的是为了解决传统现场测量手段在河流平面指标提取过程中存在的高成本低效率问题，并在保证了所提取的河道的连续性与指标数据的多样性前提下，提出了一种基于复杂网络分析技术的河流平面几何指标与拓扑结构的自动化提取方法。

本发明的技术方案为：

一种河流平面几何指标与拓扑结构的提取方法，包括以下步骤：

S1、确定输入图像，并率定初始参数。包括：

S11、输入的图像格式为卫星图片，首先确定所输入图像的分辨率，以及图像尺寸M×N，其中，M，N分别为图像的像素行数与列数；其次，按照改进的归一化差异水体指数公式：计算并储存每个图像的像素点所对应的MNDWI，其中Green与MIR分别为卫星图像中绿色波段与中红外波段的反射率。

S12、率定模型模型初始参数σ₁，σ₂与ΔW，其中σ₁与σ₂分别为去除噪声过程中膨胀与腐蚀算法所使用的卷积核的最大尺寸，而ΔW为河宽提取指标参数；初始参数的最优值均通过对所输出的河道数目的敏感性分析而确定，具体过程为，分别从0开始逐渐增加各参数值，当所输出到的河道数目趋于稳定值后，所对应的参数值即为最优值。

S2、提取河道水体并消除噪声。包括：

S21、提取河道水体：根据S11计算所得到的MNDWI指数，对图像所有像素点运用OTSU阈值分割方法，找出水体像素与陆地像素的阈值T，对于MNDWI指数大于T的像素定义为水体像素，剩下非水体像素则被当作背景去除，并输出河道水体的二值图像，值为1的像素点为河道水体像素，值为0的像素点为背景像素。

S22、噪声消除：通过S21所得到的二值图像中的像素点与卷积核进行卷积运算，其中卷积核的边缘像素值设为1，中间像素值为0；首先，遍历所有水体像素点，根据膨胀算法使其分别与尺寸为n₁(n₁＝1，3，5.....σ₁)的卷积核进行卷积，若结果为0，则标记为零星水域噪声并予以消除；其次，同样地，遍历所有背景像素点，根据腐蚀算法分别使其与尺寸为n₂(n₂＝1，3，5.....σ₂)的卷积核进行卷积，若结果为4σ₂-4，则标记为对河流平面形态无主要影响的较小的江心滩或沙坝并予以填充，最终得到整个河流系统的掩膜。

S3、刻画河道中心线骨架；包括：

S31、对于步骤S22得到的河流掩膜进行Zhang-Suen快速并行细化算法，得到河流骨架，记为矩阵S＝{X_s,Y_s}，其中X_s与Y_s分别为第s个骨架像素点的行、列坐标；

S32、根据步骤S31所得到的河流骨架，遍历所有骨架像素点，分别采用消除模版与其进行卷积，所得结果若属于集合{65，5，20，80，13，22，52，133，141，54}，则删除该骨架像素点，以保证所得河流骨架为单分辨率宽度。

S4、检测河流的分汊点与端点；包括：

S41、记骨架像素点的8个邻域中同样为骨架的像素点个数为N_s，遍历所有河流骨架像素点，并检测其N_s的取值，如果N_s＝3，则判定该骨架像素点为河流分汊点；如果N_s＝1，则判定其为河流端点。

S42、记录所有得到的河流分汊点与端点坐标，并按照从河流上游到下游的顺序对这些点进行排序。

S5、构建河流的拓扑结构，提取河流的平面几何指标。包括：

S51、根据步骤S42所得到的河流分汊点与端点进行追踪，首先，按顺序从河流的一个分汊点或者端点i开始，逐个追踪并记录与其相连的骨架像素点及其坐标，直至另一个河流分汊点或端点j结束；然后，按照相同的方式追踪另外一条与i相连的河流骨架；最后，直到所有与i相连的河流骨架均被追踪完成。

S52、河流拓扑结构的构建。河流拓扑结构的数学表达通过其河网连通矩阵A_n,n＝{a_i，j}来实现，其中n为河流分汊点与端点的数目，a_i，j表示河流中第i个与第j个分汊点或端点的连通关系，a_i,j＝1表示两分汊点或端点连通，a_i,j＝0则表示不连通。首先，假设连通矩阵A_n,n为零矩阵，设在骨架追踪过程中，若检测到一条河流骨架的两端分别为点i与j并且i与j不同时为河流端点，则a_i,j＝1,否则a_i,j＝0，最后，直至所有的河流分汊点与端点均已追踪完成，则所输出的连通矩阵即为该河流的网络拓扑结构表示。

S53、河流平面几何指标的提取。在从分汊点或端点i到j的河流骨架追踪的过程中，通过对于骨架元素点坐标运算来提取河流的河道长度与曲率，根据表达式与分别求出该骨架所代表的河道的曲线长度/>直线长度/>以及曲率/>其中/>与/>分别为连接点i和j的河流骨架中第k个像素点的横纵坐标；而X_i/Y_i与X_j/Y_j分别是河流分汊点或端点i与j的横纵坐标；K为骨架中所追踪到的像素点总个数。

求出河道的宽度：首先，在追踪过程中，根据表达式来计算连接点i与j的河流骨架中的第k个像素点的局部坡度，垂直于该坡度在河流掩膜中向该像素点左右两侧做垂线得到左右横截面，若其左右横截面差值的绝对值小于ΔW，则该横截面被判定为有效横截面，其长度被记录为像素点k的局部河宽/>然后，当一条河流骨架被追踪完成后，计算所有被记录的该骨架上的有效横截面宽度的平均值，即为该骨架所对应的河道的河宽。

本发明的有益效果是：本发明包括河道水体的提取，河流中心线骨架的刻画，河流分汊点与端点的识别与检测，及河流网络拓扑结构与平面几何指标的提取四个主要部分，节约了传统现场测量手段所带来的大量的人力物力消耗，不仅保证了所提取的河道的连续与指标的多样性，同时还考虑了江心滩，沙坝等对于河道平面形态所产生的影响，因此在辫状河流，畸形河湾等具有复杂平面形态的条件下同样具有较强的实用性。

附图说明

图1为本发明提供的一种河流平面几何指标与拓扑结构的提取方法流程图；

图2为本发明所提取的河流水体及膨胀与腐蚀算法的卷积核示意图；

图3为本发明细化算法消除模版示意图；

图4为本发明河流骨架与分汊点端点示意图；

图5为本发明河流平面几何指标提取追踪程序示意图。

具体实施方式

下面结合附图对本发明的实施例作进一步的说明。

本发明提供了一种河流平面几何指标与拓扑结构的提取方法，由河道水体的提取、河流中心线骨架的刻画，河流分汊点与端点的识别与检测和河流网络拓扑结构与平面几何指标的提取四个部分组成，如图1所示，具体包括以下步骤：

S1、确定输入图像，并率定初始参数。

该步骤具体包括以下分步骤：

S11、输入的图像格式为卫星图片(以Landasat卫星图片为例)，首先确定所输入图像I的分辨率，以及图像尺寸M×N(其中，M，N分别为图像的像素行数与列数)；其次，按照改进的归一化差异水体指数(Modified Normalized Difference Water Index,MNDWI)公式：计算并储存每个图像的像素点所对应的MNDWI，其中Green与MIR分别为卫星图像中绿色波段与中红外波段的反射率。

S12、率定模型模型初始参数σ₁，σ₂与ΔW。其中σ₁与σ₂分别为去除噪声过程中膨胀与腐蚀算法所使用的卷积核(如图2中2所示)的最大尺寸，而ΔW为河宽提取过程相关参数。初始参数的最优值均可通过对所输出的河道数目的敏感性分析而确定，具体过程为，分别从0开始逐渐增加各参数值(例如σ₁＝0，1，2......)，当所输出到的河道数目趋于稳定值后，所对应的参数值即为最优值。另外，也可采用σ₁，σ₂与ΔW的参考值，分别为图像尺寸M和N中的较小值，3倍的分辨率宽度，5倍的分辨率宽度。

S2、提取河道水体并消除噪声。

该步骤具体包括以下分步骤：

S21、提取河道水体。根据S11计算所得到的MNDWI指数，对图像所有像素点运用OTSU阈值分割方法，找出水体像素与陆地像素的阈值T，即对于MNDWI指数大于T的像素定义为水体像素，剩下非水体像素则被当作背景去除，并输出河道水体的二值图像(值为1的像素点为河道水体像素，值为0的像素点为背景像素)。

S22、噪声消除。噪声消除是通过S21所得到的二值图像中的像素点与卷积核进行卷积运算来完成的，其中卷积核的边缘像素值设为1，中间像素值为0(如图2中2所示)。首先，遍历所有水体像素点，根据膨胀算法使其分别与尺寸为n₁(n₁＝1，3，5.....σ₁)的卷积核进行卷积，若结果为0，则标记为零星水域噪声并予以消除；其次，同样地，遍历所有背景像素点，根据腐蚀算法分别使其与尺寸为n₂(n₂＝1，3，5.....σ₂)的卷积核进行卷积，若结果为(4σ₂-4)，则标记为对河流平面形态无主要影响的较小的江心滩或沙坝并予以填充，最终得到了整个河流系统的掩膜(如图2中1所示)。

S3、刻画河道中心线骨架。

该步骤具体包括以下分步骤：

S31、对于步骤S22得到的河流掩膜进行Zhang-Suen快速并行细化算法，得到河流骨架，记为矩阵S＝{X_s,Y_s}，其中X_s与Y_s分别为第s个骨架像素点的行，列坐标。

S32、根据步骤S31所得到的河流骨架，遍历所有骨架像素点，分别采用如图3所示的10种消除模版与其进行卷积，所得结果若属于集合{65，5，20，80，13，22，52，133，141，54}，则删除该骨架像素点，以保证所得河流骨架为单分辨率宽度。

S4、检测河流的分汊点与端点。

该步骤具体包括以下分步骤：

S41、记骨架像素点(如图4中3所示)的8个邻域中同样为骨架的像素点个数为N_s，遍历所有河流骨架像素点，并检测其N_s的取值，如果N_s＝3，则判定该骨架像素点为河流分汊点(如图4中4所示)；如果N_s＝1，则判定其为河流端点(如图4中5所示)。

S42、记录所有得到的河流分汊点与端点坐标，并按照从河流上游到下游的顺序对这些点进行排序。

S5、构建河流的拓扑结构，提取河流的平面几何指标。

该步骤具体包括以下分步骤：

S51、根据步骤S42所得到的河流分汊点与端点执行追踪程序。首先，该追踪程序按顺序从河流的一个分汊点或者端点i(以图5中7为例)开始，逐个追踪并记录与其相连的骨架像素点及其坐标，直至另一个河流分汊点或端点j(以图5中8为例)结束；然后，按照相同的方式追踪另外一条与i相连的河流骨架；最后，直到所有与i相连的河流骨架均被追踪完成。

S52、河流拓扑结构的构建。河流拓扑结构的数学表达通过其河网连通矩阵A_n,n＝{a_i，j}来实现，其中n为河流分汊点与端点的数目，而a_i，j则表示河流中第i个与第j个分汊点或端点的连通关系，a_i,j＝1表示两分汊点或端点连通，a_i,j＝0则表示不连通。首先，假设连通矩阵A_n,n为零矩阵，设在骨架追踪过程中，若检测到一条河流骨架的两端分别为点i与j并且i与j不同时为河流端点，则a_i,j＝1,否则a_i,j＝0，最后，直至所有的河流分汊点与端点均已追踪完成，则所输出的连通矩阵即为该河流的网络拓扑结构表示。

S53、河流平面几何指标的提取。在从分汊点或端点i到j的河流骨架追踪的过程中，通过对于骨架元素点坐标运算来提取河流的河道长度与曲率，根据表达式与分别求出该骨架所代表的河道的曲线长度/>直线长度/>以及曲率/>其中/>与/>分别为连接点i和j的河流骨架中第k个像素点的横纵坐标；而X_i/Y_i与X_j/Y_j分别是河流分汊点或端点i与j的横纵坐标；K为骨架中所追踪到的像素点总个数。

河道的宽度则是通过建立相应的横截面(如图5中6所示)所得到。类似地，首先，在追踪过程中，根据表达式来计算连接点i与j的河流骨架中的第k个像素点的局部坡度，垂直于该坡度在河流掩膜中向该像素点左右两侧做垂线得到左右横截面，若其左右横截面差值的绝对值小于ΔW，则该横截面被判定为有效横截面，其长度被记录为像素点k的局部河宽/>然后，当一条河流骨架被追踪完成后，计算所有被记录的该骨架上的有效横截面宽度的平均值，即为该骨架所对应的河道的河宽。

S52、河流平面几何指标的可视化与输出。类似于河流拓扑结构的连通矩阵，河流的平面几何指标，包括河道长度，宽度，曲率分别通过可视化矩阵L_n,n＝{l_i，j}，W_n,n＝{w_i，j}与S_n,n＝{s_i，j}输出，其中l_i，j，w_i，j与s_i，j分别表示连接河流第i个与第j个分汊点或端点的河道的河长，河宽与曲率。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (2)

1.一种河流平面几何指标与拓扑结构的提取方法，其特征在于，包括以下步骤：

S1、确定输入图像，并率定初始参数，包括：

S11、输入的图像格式为卫星图片，首先确定所输入图像的分辨率，以及图像尺寸M×N，其中，M，N分别为图像的像素行数与列数；其次，按照改进的归一化差异水体指数公式：计算并储存每个图像的像素点所对应的MNDWI，其中Green与MIR分别为卫星图像中绿色波段与中红外波段的反射率；

S12、率定模型初始参数σ₁，σ₂与ΔW，其中σ₁与σ₂分别为去除噪声过程中膨胀与腐蚀算法所使用的卷积核的最大尺寸，而ΔW为河宽提取指标参数；初始参数的最优值均通过对所输出的河道数目的敏感性分析而确定，具体过程为，分别从0开始逐渐增加各参数值，当所输出到的河道数目趋于稳定值后，所对应的参数值即为最优值；

S2、提取河道水体并消除噪声，包括：

S21、提取河道水体：根据S11计算所得到的MNDWI指数，对图像所有像素点运用OTSU阈值分割方法，找出水体像素与陆地像素的阈值T，对于MNDWI指数大于T的像素定义为水体像素，剩下非水体像素则被当作背景去除，并输出河道水体的二值图像，值为1的像素点为河道水体像素，值为0的像素点为背景像素；

S22、噪声消除：通过S21所得到的二值图像中的像素点与卷积核进行卷积运算，其中卷积核的边缘像素值设为1，中间像素值为0；首先，遍历所有水体像素点，根据膨胀算法使其分别与尺寸为n₁(n₁＝1，3，5.....σ₁)的卷积核进行卷积，若结果为0，则标记为零星水域噪声并予以消除；其次，同样地，遍历所有背景像素点，根据腐蚀算法分别使其与尺寸为n₂(n₂＝1，3，5.....σ₂)的卷积核进行卷积，若结果为4σ₂-4，则标记为对河流平面形态无主要影响的较小的江心滩或沙坝并予以填充，最终得到整个河流系统的掩膜；

S3、刻画河道中心线骨架；包括：

S31、对于步骤S22得到的河流掩膜进行Zhang-Suen快速并行细化算法，得到河流骨架，记为矩阵S＝{X_s,Y_s}，其中X_s与Y_s分别为第s个骨架像素点的行、列坐标；

S32、根据步骤S31所得到的河流骨架，遍历所有骨架像素点，分别采用消除模版与其进行卷积，所得结果若属于集合{65，5，20，80，13，22，52，133，141，54}，则删除该骨架像素点，以保证所得河流骨架为单分辨率宽度；

S4、检测河流的分汊点与端点；包括：

S41、记骨架像素点的8个邻域中同样为骨架的像素点个数为N_s，遍历所有河流骨架像素点，并检测其N_s的取值，如果N_s＝3，则判定该骨架像素点为河流分汊点；如果N_s＝1，则判定其为河流端点；

S42、记录所有得到的河流分汊点与端点坐标，并按照从河流上游到下游的顺序对这些点进行排序；

S5、构建河流的拓扑结构，提取河流的平面几何指标，包括：

S51、根据步骤S42所得到的河流分汊点与端点进行追踪，首先，按顺序从河流的一个分汊点或者端点i开始，逐个追踪并记录与其相连的骨架像素点及其坐标，直至另一个河流分汊点或端点j结束；然后，按照相同的方式追踪另外一条与i相连的河流骨架；最后，直到所有与i相连的河流骨架均被追踪完成；

S52、河流拓扑结构的构建，河流拓扑结构的数学表达通过其河网连通矩阵A_n,n＝{a_i，j}来实现，其中n为河流分汊点与端点的数目，a_i，j表示河流中第i个与第j个分汊点或端点的连通关系，a_i,j＝1表示两分汊点或端点连通，a_i,j＝0则表示不连通；首先，假设连通矩阵A_n,n为零矩阵，设在骨架追踪过程中，若检测到一条河流骨架的两端分别为点i与j并且i与j不同时为河流端点，则a_i,j＝1,否则a_i,j＝0，最后，直至所有的河流分汊点与端点均已追踪完成，则所输出的连通矩阵即为该河流的网络拓扑结构表示；

S53、河流平面几何指标的提取，在从分汊点或端点i到j的河流骨架追踪的过程中，通过对于骨架元素点坐标运算来提取河流的河道长度与曲率，根据表达式与分别求出该骨架所代表的河道的曲线长度/>直线长度/>以及曲率/>其中/>与/>分别为连接点i和j的河流骨架中第k个像素点的横纵坐标；而X_i/Y_i与X_j/Y_j分别是河流分汊点或端点i与j的横纵坐标；K为骨架中所追踪到的像素点总个数；

求出河道的宽度：首先，在追踪过程中，根据表达式来计算连接点i与j的河流骨架中的第k个像素点的局部坡度，垂直于该坡度在河流掩膜中向该像素点左右两侧做垂线得到左右横截面，若其左右横截面差值的绝对值小于ΔW，则该横截面被判定为有效横截面，其长度被记录为像素点k的局部河宽/>然后，当一条河流骨架被追踪完成后，计算所有被记录的该骨架上的有效横截面宽度的平均值，即为该骨架所对应的河道的河宽。

2.根据权利要求1所述的一种河流平面几何指标与拓扑结构的提取方法，其特征在于，还包括以下步骤：

S5、构建河流的拓扑结构，提取河流的平面几何指标。