CN110688961B - Method and system for extracting topology information of river network - Google Patents

Method and system for extracting topology information of river network Download PDF

Info

Publication number
CN110688961B
CN110688961B CN201910937320.2A CN201910937320A CN110688961B CN 110688961 B CN110688961 B CN 110688961B CN 201910937320 A CN201910937320 A CN 201910937320A CN 110688961 B CN110688961 B CN 110688961B
Authority
CN
China
Prior art keywords
water body
river
pixel
condition
image
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910937320.2A
Other languages
Chinese (zh)
Other versions
CN110688961A (en
Inventor
陈曦
崔要奎
万玮
雷添杰
李欢
沈心一
李京
方唯振
刘宝剑
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Peking University
Original Assignee
Peking University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Peking University filed Critical Peking University
Priority to CN201910937320.2A priority Critical patent/CN110688961B/en
Publication of CN110688961A publication Critical patent/CN110688961A/en
Application granted granted Critical
Publication of CN110688961B publication Critical patent/CN110688961B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/182—Network patterns, e.g. roads or rivers
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2411—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on the proximity to a decision surface, e.g. support vector machines
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/243—Classification techniques relating to the number of classes
    • G06F18/24323—Tree-organised classifiers
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/50—Depth or shape recovery
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Evolutionary Biology (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Pure & Applied Mathematics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Mathematical Optimization (AREA)
  • Computing Systems (AREA)
  • Multimedia (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Image Processing (AREA)

Abstract

The invention discloses a method and a system for extracting topology information of a river network. The method comprises the following steps: extracting the water body in the remote sensing image to obtain a water body image; extracting a river skeleton in the water body image by adopting a morphological thinning algorithm to obtain a river skeleton image; determining end points and branch points in the river skeleton image; respectively taking each endpoint as a starting point, and searching adjacent branch points of each endpoint along the skeleton path by adopting a path search algorithm to obtain a search path and a search path length corresponding to the endpoint; calculating the nearest Euclidean distance of the offshore of the adjacent nodes; deleting the search path meeting the pruning condition in the river skeleton image to obtain a pruned river skeleton image; the pruning condition is determined by the length of a search path corresponding to an endpoint, the nearest Euclidean distance of the offshore, a preset proportional coefficient threshold value and a preset absolute length threshold value; and determining a connectivity matrix of the pruned river skeleton image. The invention can realize automatic, complete and accurate extraction of the topology information of the river network.

Description

Method and system for extracting topology information of river network
Technical Field
The invention relates to the technical field of water resource monitoring and information extraction, in particular to a method and a system for extracting topology information of a river network.
Background
The surface water sources of China are various, rivers and lakes are numerous, the north-south difference is large, and the annual and annual changes are large. The identification and extraction of water body information based on a remote sensing technology become a main method for realizing large-range and high-efficiency water body mapping. After the water coverage area mapping is completed, river network topology information is further extracted, and basic data guarantee is provided for applications such as river and lake resource dynamic monitoring, water resource management right confirming and hydrological simulation.
The river network topology information mainly refers to the communication relationship between the river channels. If a curve segment is used for abstract representation of a section of non-branched river channel, the topological relation of the river segment can be described through the connection relation between the curve segments. Further, if each curve segment representing the river segment is abstracted into one element, the connectivity matrix can be used to represent the connection relationship between the curve segments, i.e. to represent the topology information of the river channel.
At present, the method for extracting the topology structure of the surface river network is mainly an artificial drawing method, for example, a surface river network database based on crowd-sourced data provided by OpenStreetMap is extracted and uploaded by a large number of volunteers in different ways (including positioning by using an on-site GPS, manual drawing based on remote sensing images, and the like), although the precision is high, some covered fine water channels can be marked, the drawing capability requirement on the volunteers providing data is high, the data quality cannot be guaranteed, data loss in some remote areas is large, the data updating is slow, and the description capability on the situations of river seasonal diversion and the like is weak.
For the automatic extraction of the topology structure of the surface river network, the existing automatic extraction method has certain limitations, and particularly has poor holding capacity for the topology structure of the river network. Several methods or software for automatically extracting river network structure information are published internationally at present, including: RivWidth software, RivaMap tool, PyRIS tool, GeoNet tool and RivWidthCloud tool. The above-mentioned tool is mainly used for extracting non-topological information of the river channel, for example: river width, river centerline, river direction of large rivers, etc. For the extraction of the river channel topological information, the RivWidth software can not only distort the extracted river channel center line, but also can not completely maintain the topological property of the river network; although the RivaMap tool can keep the central lines of most river channels, the central lines of the intersected river channels are not connected, so that the topological property of the river network cannot be completely kept; the PyRIS tool also has the defects of distorted river channel center line extraction and incapability of completely keeping the topological property of the river network; the GeoNet tool extracts the river channel based on the surface elevation data, and the extraction result of the GeoNet tool is different from the real river channel structure in the remote sensing image to a certain extent, so that the defect of extraction distortion of the center line of the river channel exists; the river channel central line extracted by the RivWidthCloud tool has a great amount of non-communication situations, the details of the central line of a part of branch river channels are lost, and the topological property of the river network cannot be completely maintained. In summary, the existing method for extracting the topology information of the river network simplifies the structure of the river network, so that the center line of the river channel is distorted (possibly passing through a land area), and the accuracy is poor, or the generated river channel has poor connectivity, so that a large number of unconnected river channel center lines are formed, and the integrity of the topology structure of the river network cannot be maintained.
Disclosure of Invention
Therefore, there is a need to provide a method and a system for extracting river network topology information, so as to realize complete and accurate extraction of the river network topology information.
In order to achieve the purpose, the invention provides the following scheme:
a method for extracting topology information of a river network comprises the following steps:
extracting the water body in the remote sensing image to obtain a water body image; the water body image consists of a plurality of water body pixel points and non-water body pixel points;
extracting a river skeleton in the water body image by adopting a morphological thinning algorithm to obtain a river skeleton image; the river skeleton image is composed of a plurality of skeleton pixel points;
determining end points and branch points in the current river skeleton image under the current iteration number m; only one pixel point in the neighborhood of eight pixels around the end point is a skeleton pixel point; more than two pixel points in the neighborhood of eight pixels around the branch point are skeleton pixel points;
respectively taking each end point as a starting point, and searching adjacent nodes of each end point along a skeleton path by adopting a path search algorithm to obtain a search path and a search path length corresponding to the end point; the adjacent node is a first node which is searched along the skeleton pixel point from the end point; the nodes are branch points or end points; the search path is composed of all skeleton pixel points passing between the end point and the adjacent node in the search process; the length of the search path is the sum of distances between all adjacent skeleton pixel points on the search path;
calculating a nearest Euclidean distance offshore of the neighboring node; the nearest Euclidean distance of the offshore is the shortest Euclidean distance between the adjacent nodes and the non-water body pixel points;
judging whether a search path meeting the pruning condition exists in the current river skeleton image under the current iteration number m;
if yes, deleting the search path meeting the pruning condition in the current river skeleton image to obtain the river skeleton image after the path deletion, adding 1 to the current iteration number m, taking the river skeleton image after the path deletion as the current river skeleton image, and returning to the determined current iteration number m to determine the end point and the branch point in the current river skeleton image; the pruning condition is determined by the search path length corresponding to the end point, the nearest Euclidean distance of the offshore, a preset proportionality coefficient threshold value and a preset absolute length threshold value; the preset proportionality coefficient threshold is a set ratio of the search path length to the nearest Euclidean distance offshore; the preset absolute length threshold is a set absolute length value of the search path;
if not, determining the river skeleton image obtained by deleting the path under the last iteration number m-1 as the pruned river skeleton image;
determining a connectivity matrix of the pruned river skeleton image; the total line number and the total column number of the connectivity matrix are equal to the sum of the number of end points and branch points in the pruned river skeleton image; the connectivity matrix represents the extracted river network topology information.
Optionally, the method for extracting the river skeleton in the water body image by using the morphological thinning algorithm to obtain the river skeleton image specifically includes:
under the current iteration number n, judging whether water body pixel points meeting a first constraint condition or a second constraint condition exist in the current water body image by adopting a morphological thinning algorithm; the first constraint condition and the second constraint condition are both determined by pixel values in an eight-pixel neighborhood around a water body pixel point;
if so, deleting the water body pixel points meeting the first constraint condition or the second constraint condition, updating the water body image, adding 1 to the current iteration number n, taking the updated water body image as the current water body image, and returning to the current iteration number n to judge whether the water body pixel points meeting the first constraint condition or the second constraint condition exist in the current water body image by adopting a morphological thinning algorithm;
and if not, taking the water body image updated under the last iteration number n-1 as a river skeleton image.
Optionally, the first constraint includes a first constraint sub-condition, a second constraint sub-condition and a third constraint sub-condition; the second constraint includes a first constraint sub-condition, a second constraint sub-condition, and a fourth constraint sub-condition;
the first constraint sub-condition is XH(p) ═ 1; wherein the content of the first and second substances,
Figure BDA0002221923700000041
p represents a water body pixel point; biA neighborhood pixel operator representing pixel p;
Figure BDA0002221923700000042
wherein x is2i-1、x2iAnd x2i+1The method comprises the steps of representing neighborhood pixel values in eight-pixel neighborhood around a pixel p, wherein i is 1,2,3 and 4, 8 pixel points in the neighborhood of a water body pixel point p are arranged in sequence anticlockwise by taking an adjacent water body pixel point on the right side of p as a starting point, and the pixel value of the adjacent pixel point on the right side of p is x1Sequentially arranged in a counterclockwise directionThe pixel values of the adjacent water body pixel points are x respectively2，x3......x8，x9Wherein x is9＝x1(ii) a The pixel value of the water body pixel point is represented by a numerical value 1, and the pixel value of the non-water body pixel point is represented by a numerical value 0;
the second constraint sub-condition is that min { n is more than or equal to 21(p),n2(p) } is less than or equal to 3; wherein the content of the first and second substances,
Figure BDA0002221923700000043
v represents a logical OR operation, and when k is 4, x2k+1＝x1；
The third constraint sub-condition is
Figure BDA0002221923700000044
Λ represents a logical and operation,
Figure BDA0002221923700000045
representing a pixel value x8The negation value of (a);
the fourth constraint sub-condition is
Figure BDA0002221923700000046
Figure BDA0002221923700000047
Representing a pixel value x4The inverted value of (d).
Optionally, the pruning condition includes a first pruning constraint condition and a second pruning constraint condition;
the first pruning constraint condition is DR(h,j)＜r·DB(j) (ii) a Wherein D isR(h, j) represents the length of the search path corresponding to the endpoint h, and j represents the adjacent node of the endpoint h; r represents a preset scaling factor threshold; dB(j) Represents the nearest euclidean distance offshore of the neighboring node j;
the second pruning constraint condition is DR(h, j) < d; where d represents a preset absolute length threshold.
Optionally, the determining a connectivity matrix of the pruned river skeleton image specifically includes:
determining an initial connectivity matrix of the pruned river skeleton image; the initial connectivity matrix is a w-row and w-column all-zero matrix, and w is the total number of nodes in the pruned river skeleton image; the nodes are end points or branch points;
for each node e in the pruned river skeleton image, searching a node f adjacent to the node e in any direction in the neighborhood of the node e by adopting a region growing method; e and f are the node numbers in the pruned river skeleton image; e is more than or equal to 1 and less than or equal to w, f is more than or equal to 1 and less than or equal to w, and e is not equal to f;
judging whether a path connecting the node e and the node f meets a path connection condition or not; the path connection condition is that all the pixels on the path are skeleton pixels, and nodes except the node e and the node f do not exist on the path;
and if so, setting the element of the e-th row and the f-th column in the initial connectivity matrix as 1.
The invention also provides a river network topology information extraction system, which comprises:
the water body image extraction module is used for extracting the water body in the remote sensing image to obtain a water body image; the water body image consists of a plurality of water body pixel points and non-water body pixel points;
the skeleton image extraction module is used for extracting a river skeleton in the water body image by adopting a morphological thinning algorithm to obtain a river skeleton image; the river skeleton image is composed of a plurality of skeleton pixel points;
the first determining module is used for determining end points and branch points in the current river skeleton image under the current iteration number m; only one pixel point in the neighborhood of eight pixels around the end point is a skeleton pixel point; more than two pixel points in the neighborhood of eight pixels around the branch point are skeleton pixel points;
the path searching module is used for respectively taking each end point as a starting point and searching adjacent nodes of each end point along a skeleton path by adopting a path searching algorithm to obtain a searching path and a searching path length corresponding to the end point; the adjacent node is a first node which is searched along the skeleton pixel point from the end point; the nodes are branch points or end points; the search path is composed of all skeleton pixel points passing between the end point and the adjacent node in the search process; the length of the search path is the sum of distances between all adjacent skeleton pixel points on the search path;
an offshore euclidean distance calculation module for calculating the nearest offshore euclidean distance of the neighboring nodes; the nearest Euclidean distance of the offshore is the shortest Euclidean distance between the adjacent nodes and the non-water body pixel points;
the judging module is used for judging whether a search path meeting the pruning condition exists in the current river skeleton image under the current iteration number m;
the pruning module is used for deleting the search path meeting the pruning condition in the current river skeleton image if the search path meeting the pruning condition exists in the current river skeleton image under the previous iteration number m to obtain the river skeleton image after the path deletion, adding 1 to the current iteration number m, taking the river skeleton image after the path deletion as the current river skeleton image, and returning to the first determining module; the pruning condition is determined by the search path length corresponding to the end point, the nearest Euclidean distance of the offshore, a preset proportionality coefficient threshold value and a preset absolute length threshold value; the preset proportionality coefficient threshold is a set ratio of the search path length to the nearest Euclidean distance offshore; the preset absolute length threshold is a set absolute length value of the search path;
the second determining module is used for determining the river skeleton image obtained by deleting the path under the previous iteration number m-1 as the pruned river skeleton image if the current river skeleton image does not have the search path meeting the pruning condition under the previous iteration number m;
the third determination module is used for determining a connectivity matrix of the pruned river skeleton image; the total line number and the total column number of the connectivity matrix are equal to the sum of the number of end points and branch points in the pruned river skeleton image; the connectivity matrix represents the extracted river network topology information.
Optionally, the skeleton image extraction module specifically includes:
the first judgment unit is used for judging whether water body pixel points meeting a first constraint condition or a second constraint condition exist in the current water body image by adopting a morphological thinning algorithm under the current iteration number n; the first constraint condition and the second constraint condition are both determined by pixel values in an eight-pixel neighborhood around a water body pixel point;
the water body image updating unit is used for deleting the water body pixel points meeting the first constraint condition or the second constraint condition if the water body pixel points meeting the first constraint condition or the second constraint condition exist in the current water body image, updating the water body image, adding 1 to the current iteration number n, taking the updated water body image as the current water body image, and returning the current water body image to the first judging unit;
and the skeleton image determining unit is used for taking the water body image updated under the last iteration number n-1 as a river skeleton image if no water body pixel points meeting the first constraint condition or the second constraint condition exist in the current water body image.
Optionally, the first constraint includes a first constraint sub-condition, a second constraint sub-condition and a third constraint sub-condition; the second constraint includes a first constraint sub-condition, a second constraint sub-condition, and a fourth constraint sub-condition;
the first constraint sub-condition is XH(p) ═ 1; wherein the content of the first and second substances,
Figure BDA0002221923700000061
p represents a water body pixel point; biA neighborhood pixel operator representing pixel p;
Figure BDA0002221923700000062
wherein x is2i-1、x2iAnd x2i+1The method comprises the steps of representing neighborhood pixel values in eight-pixel neighborhood around a pixel p, wherein i is 1,2,3 and 4, 8 pixel points in the neighborhood of a water body pixel point p are arranged in sequence anticlockwise by taking an adjacent water body pixel point on the right side of p as a starting point, and the pixel value of the adjacent pixel point on the right side of p is x1The pixel values of the adjacent water body pixel points which are sequentially arranged anticlockwise are x respectively2，x3......x8，x9Wherein x is9＝x1(ii) a The pixel value of the water body pixel point is represented by a numerical value 1, and the pixel value of the non-water body pixel point is represented by a numerical value 0;
the second constraint sub-condition is that min { n is more than or equal to 21(p),n2(p) } is less than or equal to 3; wherein the content of the first and second substances,
Figure BDA0002221923700000071
v represents a logical OR operation, and when k is 4, x2k+1＝x1；
The third constraint sub-condition is
Figure BDA0002221923700000072
Λ represents a logical and operation,
Figure BDA0002221923700000073
representing a pixel value x8The negation value of (a);
the fourth constraint sub-condition is
Figure BDA0002221923700000074
Figure BDA0002221923700000075
Representing a pixel value x4The inverted value of (d).
Optionally, the pruning condition includes a first pruning constraint condition and a second pruning constraint condition;
the first pruning constraint condition is DR(h,j)＜r·DB(j) (ii) a Wherein D isR(h, j) represents the length of the search path corresponding to the endpoint h, and j represents the adjacent node of the endpoint h; r represents a preset scaling factor threshold; dB(j) Represents the nearest euclidean distance offshore of the neighboring node j;
the second pruning constraint condition is DR(h, j) < d; where d represents a preset absolute length threshold.
Optionally, the third determining module specifically includes:
the initial matrix determining unit is used for determining an initial connectivity matrix of the pruned river skeleton image; the initial connectivity matrix is a w-row and w-column all-zero matrix, and w is the total number of nodes in the pruned river skeleton image; the nodes are end points or branch points;
the searching unit is used for searching a node f adjacent to the node e in any direction in the neighborhood of the node e by adopting a region growing method for each node e in the pruned river skeleton image; e and f are the node numbers in the pruned river skeleton image; e is more than or equal to 1 and less than or equal to w, f is more than or equal to 1 and less than or equal to w, and e is not equal to f;
the second judging unit is used for judging whether a path connecting the node e and the node f meets the path connecting condition or not; the path connection condition is that all the pixels on the path are skeleton pixels, and nodes except the node e and the node f do not exist on the path;
and the communication matrix determining unit is used for setting the element of the e-th row and the f-th column in the initial communication matrix as 1 if the path connecting the node e and the node f meets the path connecting condition.
Compared with the prior art, the invention has the beneficial effects that:
the invention provides a method and a system for extracting topology information of a river network. According to the method, a morphological thinning algorithm is adopted to extract a river skeleton in a water body image to obtain a river skeleton image, and river network topology information is extracted based on a morphological method, so that the integrity of a river network topology structure is guaranteed, and the connectivity and the integrity of the river network topology are better; in the method, each endpoint is taken as a starting point, a path search algorithm is adopted to search the nearest branch point of each endpoint along a skeleton path, and calculating the nearest Euclidean distance of the nearest branch point, deleting the search path meeting the pruning condition in the river skeleton image to obtain the pruned river skeleton (river channel center line) image, wherein the pruning condition is determined by the search path length corresponding to the end point, the nearest Euclidean distance off the shore, a preset proportionality coefficient threshold value and a preset absolute length threshold value, the method for trimming the center line branches of the river channel has controllability, not only ensures that real river network branches are not lost, but also effectively reduces residual false river channel branches, the extracted river network has high structural integrity and low redundancy, and the extracted river channel center line does not pass through a land area or omit a branch river channel due to the existence of a braided river network or a river center island, so that the accuracy is high; in addition, the topology of the curve segment is expressed by a connectivity matrix in the method, and the automatic extraction of the river channel topology information is realized.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings needed to be used in the embodiments will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and it is obvious for those skilled in the art to obtain other drawings without inventive exercise.
Fig. 1 is a flowchart of a method for extracting topology information of a river network according to embodiment 1 of the present invention;
FIG. 2 is a schematic diagram of positions of a water pixel p and pixels in a neighborhood of eight pixels around the water pixel p in a morphological refinement algorithm according to embodiment 2 of the present invention;
FIG. 3 is a schematic diagram of branch points and end points of the skeleton according to embodiment 2 of the present invention;
FIG. 4 is a schematic diagram of the skeleton pruning process in example 2 of the present invention;
fig. 5 is a graph showing the pruning results for r-0.8 and r-1.2 in example 2 of the present invention;
fig. 6 is a graph showing the pruning results for r ═ 1.8 and r ═ 2.0 in example 2 of the present invention;
fig. 7 is a graph showing the pruning results for r-3.0 and r-5.0 in example 2 of the present invention;
fig. 8 is a schematic structural diagram of a river topology information extraction system according to embodiment 3 of the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
In order to make the aforementioned objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in further detail below.
Fig. 1 is a flowchart of a method for extracting topology information of a river network according to embodiment 1 of the present invention. Referring to fig. 1, the method for extracting topology information of a river network according to this embodiment includes:
step S1: and extracting the water body in the remote sensing image to obtain a water body image.
The water body image is composed of a plurality of water body pixel points and non-water body pixel points, wherein the pixel values of the water body pixel points are represented by a numerical value 1, and the pixel values of the non-water body pixel points are represented by a numerical value 0.
Step S2: and extracting the river skeleton in the water body image by adopting a morphological thinning algorithm to obtain a river skeleton image. The river skeleton image is composed of a plurality of skeleton pixel points.
The step S2 specifically includes:
under the current iteration number n, judging whether water body pixel points meeting a first constraint condition or a second constraint condition exist in the current water body image by adopting a morphological thinning algorithm; the first constraint condition and the second constraint condition are both determined by pixel values in an eight-pixel neighborhood around a water body pixel point.
And if so, deleting the water body pixel points meeting the first constraint condition or the second constraint condition, updating the water body image, adding 1 to the current iteration number n, taking the updated water body image as the current water body image, and returning to the current iteration number n to judge whether the water body pixel points meeting the first constraint condition or the second constraint condition exist in the current water body image by adopting a morphological thinning algorithm.
And if not, taking the water body image updated under the last iteration number n-1 as a river skeleton image.
As an optional implementation manner, the first constraint condition includes a first constraint sub-condition G1, a second constraint sub-condition G2, and a third constraint sub-condition G3, and the water body pixel points satisfying the first constraint condition are the water body pixel points satisfying the first constraint sub-condition G1, the second constraint sub-condition G2, and the third constraint sub-condition G3 at the same time; the second constraint condition comprises a first constraint sub-condition G1, a second constraint sub-condition G2 and a fourth constraint sub-condition G4, and the water body pixel points meeting the second constraint condition are the water body pixel points meeting the first constraint sub-condition G1, the second constraint sub-condition G2 and the fourth constraint sub-condition G4 at the same time.
Specifically, the first constraint sub-condition G1: xH(p) ═ 1; wherein the content of the first and second substances,
Figure BDA0002221923700000101
p represents a water body pixel point; biA neighborhood pixel operator representing pixel p;
Figure BDA0002221923700000102
wherein x is2i-1、x2iAnd x2i+1The method comprises the steps of representing neighborhood pixel values in eight-pixel neighborhood around a pixel p, wherein i is 1,2,3 and 4, 8 pixel points in the neighborhood of a water body pixel point p are arranged in sequence anticlockwise by taking an adjacent water body pixel point on the right side of p as a starting point, and the pixel value of the adjacent pixel point on the right side of p is x1The pixel values of the adjacent water body pixel points which are sequentially arranged anticlockwise are x respectively2，x3......x8，x9Wherein x is9＝x1。
Second constraint sub-condition G2: min { n is not less than 21(p),n2(p) } is less than or equal to 3; wherein the content of the first and second substances,
Figure BDA0002221923700000103
Figure BDA0002221923700000104
v represents a logical OR operation, and when k is 4, x2k+1＝x1。
Third constraint sub-condition G3:
Figure BDA0002221923700000105
Λ represents a logical and operation,
Figure BDA0002221923700000106
representing a pixel value x8The inverted value of (d).
Fourth constraint sub-condition G4:
Figure BDA0002221923700000107
Figure BDA0002221923700000108
representing a pixel value x4The inverted value of (d).
Step S3: and determining the end points and branch points in the current river skeleton image under the current iteration number m.
Only one pixel point in the neighborhood of eight pixels around the end point is a skeleton pixel point; and more than two pixel points in the neighborhood of eight pixels around the branch point are skeleton pixel points.
Step S4: and respectively taking each end point as a starting point, and searching adjacent nodes of each end point along a skeleton path by adopting a path search algorithm to obtain a search path and a search path length corresponding to the end point.
The adjacent node is a first node which is searched along the skeleton pixel point from the end point; the nodes are branch points or end points; the search path is composed of all skeleton pixel points passing between the end point and the adjacent node in the search process; the length of the search path is the sum of distances between all adjacent skeleton pixel points on the search path.
Step S5: calculating a nearest Euclidean distance offshore of the neighboring node.
The nearest Euclidean distance of the offshore is the shortest Euclidean distance between the adjacent nodes and the non-water body pixel points.
Step S6: and judging whether a search path meeting the pruning condition exists in the current river skeleton image under the current iteration number m. If so, go to step S7, otherwise, go to step S8.
Step S7: deleting the search path meeting the pruning condition in the current river skeleton image to obtain the river skeleton image after the path deletion, adding 1 to the current iteration number m, taking the river skeleton image after the path deletion as the current river skeleton image, and returning to the step S3.
The pruning condition is determined by the search path length corresponding to the end point, the nearest Euclidean distance of the offshore, a preset proportionality coefficient threshold value and a preset absolute length threshold value; the preset proportionality coefficient threshold is a set ratio of the search path length to the nearest Euclidean distance offshore; the preset absolute length threshold is a set absolute length value of the search path.
As an optional implementation manner, the pruning condition includes a first pruning constraint and a second pruning constraint. The search path satisfying the pruning condition is the search path satisfying both the first pruning constraint condition and the second pruning constraint condition.
First pruning constraint: dR(h,j)＜r·DB(j) (ii) a Wherein D isR(h, j) represents the length of the search path corresponding to the endpoint h, and j represents the adjacent node of the endpoint h; r represents a preset scaling factor threshold; dB(j) Representing the nearest euclidean distance offshore of the neighboring node j.
Second pruning constraint conditions: dR(h, j) < d; where d represents a preset absolute length threshold.
Step S8: and determining the river skeleton image obtained by deleting the path under the last iteration number m-1 as the pruned river skeleton image.
Step S9: and determining a connectivity matrix of the pruned river skeleton image.
The total row number and the total column number of the connectivity matrix are the sum of the number of end points and branch points in the pruned river skeleton image; the connectivity matrix represents the extracted river network topology information.
The step S9 specifically includes:
determining an initial connectivity matrix of the pruned river skeleton image; the initial connectivity matrix is a w-row and w-column all-zero matrix, and w is the total number of nodes in the pruned river skeleton image; the nodes are end points or branch points.
For each node e in the pruned river skeleton image, searching a node f adjacent to the node e in any direction in the neighborhood of the node e by adopting a region growing method; e and f are the node numbers in the pruned river skeleton image; e is more than or equal to 1 and less than or equal to w, f is more than or equal to 1 and less than or equal to w, and e is not equal to f.
Judging whether a path connecting the node e and the node f meets a path connection condition or not; and the path connection condition is that all the pixels on the path are skeleton pixels, and nodes except the node e and the node f do not exist on the path.
And if so, setting the element of the e-th row and the f-th column in the initial connectivity matrix as 1. And when the element in the connectivity matrix is 1, the node corresponding to the row number is the adjacent node of the node corresponding to the column number.
According to the river network topology information extraction method, the river network topology information is extracted based on a morphological method, the integrity of a river network topology structure is guaranteed, and the connectivity and the integrity of the river network topology are better; the adopted river channel center line branch trimming method has controllability, not only ensures that real river network branches are not lost, but also effectively reduces residual false river channel branches, so that the extracted river network has high structural integrity and low redundancy, the extracted river channel center line does not pass through a land area or leave out branch river channels due to the existence of a braided river network or a river center island, and the accuracy is high; the topology of the curve segment is expressed by a connectivity matrix, so that the automatic extraction of the topology information of the river channel is realized; the method has good adaptability to remote sensing images with different resolutions, and the river network topological structure can be extracted by adopting the embodiment as long as the connected water body range can be accurately extracted from the remote sensing images obtained with different resolutions and different time.
The embodiment provides a more specific method for extracting river network topology information, which includes:
1) extracting water body
There are many mature methods for extracting water body based on remote sensing image. In this embodiment, the water body in the remote sensing image is extracted, and the obtained water body image (the non-zero value of the pixel in the image is used for representing the water body, and the zero value of the pixel is used for representing the non-water body) is used as the only input data in the subsequent process.
In this embodiment, the water body image extraction method may adopt a visual interpretation method, a water body index method, a surface coverage classification method, or the like. The visual interpretation method specifically comprises the steps of judging whether each pixel in an image represents a water body or not by a manual visual identification method; the water body index method is to construct a water body index through multi-band operation of a remote sensing image, and the commonly used water body index comprises normalized difference water body indexes NDWI, MNDWI, an automatic water body extraction index AWEI and the like; the earth surface coverage classification method mainly adopts the traditional supervised and unsupervised classification methods, extracts water as one of the categories in the earth surface coverage classification system, and commonly adopts the unsupervised classification methods such as K-means, IsoDATA and the like, and the supervised classification methods such as Support Vector Machine (SVM), decision tree (decision tree), K-neighbor (KNN), Random Forest (Random Forest) and the like. The water body image extraction method can be used for extracting the water body image from the remote sensing image independently or jointly.
2) Extraction of river skeleton by morphological refinement algorithm
The morphological refinement (fining) algorithm adopted in this embodiment refers to a kind of algorithm, and functions to convert the connected water body pixels into a "skeleton" with a width of only 1 pixel without changing the topological structure thereof. Assuming that R represents a water body pixel point set in a water body image, and S is an extracted water body skeleton. First, S is initialized, and S is made R. The definition of 8 neighborhood pixels around a pixel point p (p e S) representing a body of water is shown in FIG. 2. Defining 4 constraints, which are: a first constraint sub-condition G1, a second constraint sub-condition G2, a third constraint sub-condition G3, and a third constraint sub-condition G4.
First constraint sub-condition G1: xH(p) ═ 1; wherein the content of the first and second substances,
Figure BDA0002221923700000131
p represents a water body pixel point; biA neighborhood pixel operator representing pixel p;
Figure BDA0002221923700000132
wherein x is2i-1、x2iAnd x2i+1The method comprises the steps of representing neighborhood pixel values in eight-pixel neighborhood around a pixel p, wherein i is 1,2,3 and 4, 8 pixel points in the neighborhood of a water body pixel point p are arranged in sequence anticlockwise by taking an adjacent water body pixel point on the right side of p as a starting point, and the pixel value of the adjacent pixel point on the right side of p is x1The pixel values of the adjacent water body pixel points which are sequentially arranged anticlockwise are x respectively2，x3......x8，x9Wherein x is9＝x1。
Second constraint sub-condition G2: min { n is not less than 21(p),n2(p) } is less than or equal to 3; wherein the content of the first and second substances,
Figure BDA0002221923700000133
Figure BDA0002221923700000134
v represents a logical OR operation, and when k is 4, x2k+1＝x1。
Third constraint sub-condition G3:
Figure BDA0002221923700000135
Λ represents a logical and operation,
Figure BDA0002221923700000136
representing a pixel value x8The inverted value of (d).
Fourth constraint sub-condition G4:
Figure BDA0002221923700000137
Figure BDA0002221923700000138
representing a pixel value x4The inverted value of (d).
In the first layer sub-loop, there will be and only pixels p that satisfy the constraints G1, G2, and G3 deleted from the S set; in the second layer sub-loop, there will be and only pixels p that satisfy the constraints G1, G2, and G4 deleted from the S set; and repeating the first layer of sub-loop and the second layer of sub-loop until the image is not changed any more, wherein the finally obtained pixel set S is the extracted river skeleton image. After the process, the topological structure of the water body can be proved to be unchanged, but some skeleton branches are false branches caused by the roughness of the edge of the river channel and need to be removed in the next step.
3) Pruning the skeleton
Due to the existence of false branches, the skeleton needs to be pruned to obtain a simplified river skeleton for constructing topological information. In this step, 2 customized parameters are set as the judgment conditions for pruning, one of which is a preset scaling factor threshold, the set ratio of the search path length to the nearest euclidean distance is denoted as r, and the other is a preset absolute length threshold, that is, the set absolute length value of the search path is denoted as d, and the function of which will be described in the following part of this step.
31) First, branch points and end points of the skeleton are determined, as shown in fig. 3. Assuming that p represents a point on the water skeleton (p ∈ S), the endpoint N of the skeleton is definedendRefers to a point having only 1 pixel belonging to the skeleton in the 8 neighborhoods thereof, the branch point N of the skeletonbranchMeaning that there are only more than 2 pixels in its 8-neighborhood that belong to the skeleton, can be represented by the following formula:
Figure BDA0002221923700000141
p in FIG. 31Is an end point, p2Is a branch point, p3Is a common skeleton point.
32) Then, starting from each end point, the distance from the nearest branch point is searched along the skeleton path. The optional methods comprise a region growing method, a Dijkstra algorithm, a Floyid algorithm and other shortest path routing algorithms. The region growing method is characterized in that an end point h is used as a seed point, an adjacent skeleton point is searched, the skeleton point is used as the seed point, the next skeleton point is searched until the found skeleton point belongs to a branch point j belonging to NbranchRecord the path G in the whole search processR(h, j), and its length DR(h, j); dijkstra algorithm, Floyid algorithm and other shortest path routing algorithmsThe method is a classic network shortest path search algorithm, if a framework is regarded as a network, all end points and branch points are regarded as network nodes, the algorithm can calculate the shortest path from any end point to a branch point and the path length thereof, and a branch point j belonging to N and being closest to some end point h can be found from the shortest pathbranchThe shortest path G can also be obtainedR(h, j), and the length D thereofR(h,j)。
33) Then, the Euclidean distance (the nearest Euclidean distance off the shore) from each branch point to the non-water body pixel closest to the branch point is calculated and recorded as DB(j) The calculation method may be: calculating Euclidean distance according to the coordinate of a branch point j and the coordinate of a non-water body pixel on the outer edge of the water body, and selecting the minimum value as DB(j)。
34) And finally, judging whether pruning is carried out or not according to the 2 user customized parameters r and d. If and only if the path search result satisfies the following two constraints, the path G should be pruned, i.e., deleted, from the skeletonRAll pixels in (h, j) except the corresponding branch point.
Constraint 1: dR(h,j)＜r·DB(j)
Constraint 2: dR(h,j)＜d
The pruning process is shown in fig. 4, in which part (a) in fig. 4 represents a river skeleton image before pruning, part (b) in fig. 4 represents a river skeleton image after pruning, part (c) in fig. 4 represents a topology structure diagram of a river network formed after pruning, and part (d) in fig. 4 represents a topology structure diagram of a river network formed when r is 1.5RThe (h, j) path requires a partial schematic of pruning. Wherein each reference numeral in section (c) in fig. 4 represents one node (branch point or end point) in the pruned river skeleton.
Therefore, the user customized parameter r is equivalent to a proportional coefficient threshold of the length of the skeleton branch path and the width of the river channel, the branch path is reserved only if the parameter r is larger than the threshold, otherwise, the branch path is cut off, the r can be generally set to be 1.5-2.0, obvious false branches can be eliminated, and shorter branches can be eliminated if a larger value is set; the user-customized parameter d is equivalent to an absolute length threshold of the branch path length of the skeleton, and the larger the value of d is set, the longer the maximum length of the deleted branch is, in this embodiment, d is a non-negative number, and the value range of d may be from 0 to plus infinity. Generally, the larger the parameters r and d are, the more concise the retained skeleton is, the more significant the topology structure of the main path of the river is, and the influence of the values of different r on the pruning result and the calculation efficiency is as shown in fig. 5 to 7. Referring to fig. 5-7, T represents the running time under the test platform (the unit: second, the larger the value is, the longer the time is), N represents the number of main nodes reserved for the final pruning result (the larger the value is, the more the reserved nodes are), the test platform is (CPU: intel (r) core (tm) i7-7500U @2.70GHz, RAM:8.00GB, OS: Windows10 x64), the reference number in each figure represents one node (branch point or end point) in the river skeleton after pruning, and the more the reference number is, the more the reserved nodes are. Part (a) of fig. 5 shows pruning results and calculation efficiency corresponding to r-0.8, T-219.23S, N-68, and part (b) of fig. 5 shows pruning results and calculation efficiency corresponding to r-1.2, T-73.04S, N-38; part (a) of fig. 6 indicates the pruning result and the calculation efficiency corresponding to r-1.8, T-11.15S, and N-14, and part (b) of fig. 6 indicates the pruning result and the calculation efficiency corresponding to r-2.0, T-11.15S, and N-14; part (a) of fig. 7 shows the pruning result and the calculation efficiency corresponding to r-3.0, T-8.75S, and N-12, and part (b) of fig. 7 shows the pruning result and the calculation efficiency corresponding to r-5.0, T-6.51S, and N-10.
4) Constructing a river network topological relation according to the skeleton image after pruning
After the skeleton image (the river channel central line) after pruning is obtained in the step 3), in order to represent the topological relation of the river network, the communication relation of the river channel central line is required to be represented by a connectivity matrix. The connectivity matrix is a symmetric matrix C, the number of rows and columns of the symmetric matrix C is equal to the total number of nodes (the sum of the number of branch points and end points) in the pruned skeleton, and C is initially set to be an all-zero matrix.
The process can adopt a region growing method, namely, taking a node (end point or branch point) e as a starting point, gradually searching skeleton points in any direction in the neighborhood until a next node (end point or branch point) f is encountered, and recording the path of the process
Figure BDA0002221923700000161
And its length
Figure BDA0002221923700000162
If there is a corresponding value C (e, f) in the connectivity matrix equal to 0, and
Figure BDA0002221923700000163
then C (e, f) is 1. When all paths are exhausted, the construction of the connectivity matrix is completed. The connectivity matrix can represent the communication relation between the nodes in the river channel path, thereby playing a role in river network topology analysis, hydrological modeling and other applications.
According to the river network topology information extraction method, the river network topology information is extracted based on a morphological method, the integrity of a river network topology structure is guaranteed, and the connectivity and the integrity of the river network topology are better; the adopted river channel center line branch trimming method has controllability, not only ensures that real river network branches are not lost, but also effectively reduces residual false river channel branches, so that the extracted river network has high structural integrity and low redundancy, the extracted river channel center line does not pass through a land area or leave out branch river channels due to the existence of a braided river network or a river center island, and the accuracy is high; the topology of the curve segment is expressed by a connectivity matrix, so that the automatic extraction of the topology information of the river channel is realized; by setting a proper proportional coefficient threshold of the branch path length and the river channel width and an absolute length threshold of the branch path, the accuracy is further improved, and meanwhile, the calculation efficiency is also improved; the method has good adaptability to remote sensing images with different resolutions.
The invention also provides a river network topology information extraction system, and fig. 8 is a schematic structural diagram of a river network topology information extraction system in embodiment 3 of the invention. The river network topology information extraction system comprises:
the water body image extraction module 801 is used for extracting the water body in the remote sensing image to obtain a water body image; the water body image is composed of a plurality of water body pixel points and non-water body pixel points.
The skeleton image extraction module 802 is configured to extract a river skeleton in the water body image by using a morphological refinement algorithm to obtain a river skeleton image; the river skeleton image is composed of a plurality of skeleton pixel points.
A first determining module 803, configured to determine an endpoint and a branch point in the current river skeleton image under the current iteration number m; only one pixel point in the neighborhood of eight pixels around the end point is a skeleton pixel point; and more than two pixel points in the neighborhood of eight pixels around the branch point are skeleton pixel points.
A path searching module 804, configured to search, by using each endpoint as a starting point and using a path searching algorithm, a neighboring node of each endpoint along a skeleton path, so as to obtain a search path and a search path length corresponding to the endpoint; the adjacent node is a first node which is searched along the skeleton pixel point from the end point; the nodes are branch points or end points; the search path is composed of all skeleton pixel points passing between the end point and the adjacent node in the search process; the length of the search path is the sum of distances between all adjacent skeleton pixel points on the search path.
An offshore euclidean distance calculation module 805 for calculating the nearest euclidean distance offshore of the neighboring node; the nearest Euclidean distance of the offshore is the shortest Euclidean distance between the adjacent nodes and the non-water body pixel points.
The determining module 806 is configured to determine whether a search path meeting the pruning condition exists in the current river skeleton image under the current iteration number m.
A pruning module 807, configured to delete the search path satisfying the pruning condition in the current river skeleton image if the search path satisfying the pruning condition exists in the current river skeleton image for the previous iteration time m, obtain a river skeleton image after the path deletion, add 1 to the current iteration time m, use the river skeleton image after the path deletion as the current river skeleton image, and return to the first determining module 803; the pruning condition is determined by the search path length corresponding to the end point, the nearest Euclidean distance of the offshore, a preset proportionality coefficient threshold value and a preset absolute length threshold value; the preset proportionality coefficient threshold is a set ratio of the search path length to the nearest Euclidean distance offshore; the preset absolute length threshold is a set absolute length value of the search path.
The second determining module 808 is configured to determine, if the current river skeleton image does not have a search path that meets the pruning condition under the previous iteration number m, the river skeleton image from which the path obtained under the previous iteration number m-1 is deleted as the pruned river skeleton image.
A third determining module 809, configured to determine a connectivity matrix of the pruned river skeleton image; the total line number and the total column number of the connectivity matrix are equal to the sum of the number of end points and branch points in the pruned river skeleton image; the connectivity matrix represents the extracted river network topology information.
As an optional implementation manner, the skeleton image extraction module 802 specifically includes:
the first judgment unit is used for judging whether water body pixel points meeting a first constraint condition or a second constraint condition exist in the current water body image by adopting a morphological thinning algorithm under the current iteration number n; the first constraint condition and the second constraint condition are both determined by pixel values in an eight-pixel neighborhood around a water body pixel point.
And the water body image updating unit is used for deleting the water body pixel points meeting the first constraint condition or the second constraint condition if the water body pixel points meeting the first constraint condition or the second constraint condition exist in the current water body image, updating the water body image, adding 1 to the current iteration number n, taking the updated water body image as the current water body image, and returning the current water body image to the first judging unit.
And the skeleton image determining unit is used for taking the water body image updated under the last iteration number n-1 as a river skeleton image if no water body pixel points meeting the first constraint condition or the second constraint condition exist in the current water body image.
As an optional implementation, the first constraint includes a first constraint sub-condition, a second constraint sub-condition, and a third constraint sub-condition; the second constraint includes a first constraint sub-condition, a second constraint sub-condition, and a fourth constraint sub-condition.
The first constraint sub-condition is XH(p) ═ 1; wherein the content of the first and second substances,
Figure BDA0002221923700000181
p represents a water body pixel point; biA neighborhood pixel operator representing pixel p;
Figure BDA0002221923700000182
wherein x is2i-1、x2iAnd x2i+1The method comprises the steps of representing neighborhood pixel values in eight-pixel neighborhood around a pixel p, wherein i is 1,2,3 and 4, 8 pixel points in the neighborhood of a water body pixel point p are arranged in sequence anticlockwise by taking an adjacent water body pixel point on the right side of p as a starting point, and the pixel value of the adjacent pixel point on the right side of p is x1The pixel values of the adjacent water body pixel points which are sequentially arranged anticlockwise are x respectively2，x3......x8，x9Wherein x is9＝x1(ii) a The pixel value of the water body pixel point is represented by a numerical value 1, and the pixel value of the non-water body pixel point is represented by a numerical value 0;
the second constraint sub-condition is that min { n is more than or equal to 21(p),n2(p) } is less than or equal to 3; wherein the content of the first and second substances,
Figure BDA0002221923700000183
v represents a logical OR operation, and when k is 4, x2k+1＝x1；
The third constraint sub-condition is
Figure BDA0002221923700000184
Λ represents a logical and operation,
Figure BDA0002221923700000185
representing a pixel value x8The negation value of (a);
the fourth constraint sub-condition is
Figure BDA0002221923700000186
Figure BDA0002221923700000187
Representing a pixel value x4The inverted value of (d).
As an optional implementation manner, the pruning condition includes a first pruning constraint and a second pruning constraint.
The first pruning constraint condition is DR(h,j)＜r·DB(j) (ii) a Wherein D isR(h, j) represents the length of the search path corresponding to the endpoint h, and j represents the adjacent node of the endpoint h; r represents a preset scaling factor threshold; dB(j) Represents the nearest euclidean distance offshore of the neighboring node j;
the second pruning constraint condition is DR(h, j) < d; where d represents a preset absolute length threshold.
As an optional implementation manner, the third determining module specifically includes:
the initial matrix determining unit is used for determining an initial connectivity matrix of the pruned river skeleton image; the initial connectivity matrix is a w-row and w-column all-zero matrix, and w is the total number of nodes in the pruned river skeleton image; the nodes are end points or branch points.
The searching unit is used for searching a node f adjacent to the node e in any direction in the neighborhood of the node e by adopting a region growing method for each node e in the pruned river skeleton image; e and f are the node numbers in the pruned river skeleton image; e is more than or equal to 1 and less than or equal to w, f is more than or equal to 1 and less than or equal to w, and e is not equal to f.
The second judging unit is used for judging whether a path connecting the node e and the node f meets the path connecting condition or not; and the path connection condition is that all the pixels on the path are skeleton pixels, and nodes except the node e and the node f do not exist on the path.
And the communication matrix determining unit is used for setting the element of the e-th row and the f-th column in the initial communication matrix as 1 if the path connecting the node e and the node f meets the path connecting condition.
The river network topology information extraction system of the embodiment extracts the river network topology information based on a morphological method, and the extracted river network center line has good river network topology connectivity and integrity; two user-customized parameters are provided, so that the riverway center line branch trimming process is controllable, not only is the real river network branch not lost, but also the residual false riverway branches are effectively reduced, the extracted river network structure is high in integrity and low in redundancy, and logic errors such as the situation that a riverway passes through a land area can be avoided; the skeleton structure in the image is converted into a connectivity matrix structure which can be expressed by mathematics, so that the abstract extraction of the topology information of the river network is realized, the data volume is remarkably reduced, and important basic data support is provided for subsequent applications such as river network topology analysis, hydrological modeling and the like; the method has good adaptability to remote sensing images with different resolutions.
The embodiments in the present description are described in a progressive manner, each embodiment focuses on differences from other embodiments, and the same and similar parts among the embodiments are referred to each other. For the system disclosed by the embodiment, the description is relatively simple because the system corresponds to the method disclosed by the embodiment, and the relevant points can be referred to the method part for description.
The principles and embodiments of the present invention have been described herein using specific examples, which are provided only to help understand the method and the core concept of the present invention; meanwhile, for a person skilled in the art, according to the idea of the present invention, the specific embodiments and the application range may be changed. In view of the above, the present disclosure should not be construed as limiting the invention.

Claims (8)

1. A method for extracting topology information of a river network is characterized by comprising the following steps:
extracting the water body in the remote sensing image to obtain a water body image; the water body image consists of a plurality of water body pixel points and non-water body pixel points;
extracting a river skeleton in the water body image by adopting a morphological thinning algorithm to obtain a river skeleton image; the river skeleton image is composed of a plurality of skeleton pixel points;
determining end points and branch points in the current river skeleton image under the current iteration number m; only one pixel point in the neighborhood of eight pixels around the end point is a skeleton pixel point; more than two pixel points in the neighborhood of eight pixels around the branch point are skeleton pixel points;
respectively taking each end point as a starting point, and searching adjacent nodes of each end point along a skeleton path by adopting a path search algorithm to obtain a search path and a search path length corresponding to the end point; the adjacent node is a first node which is searched along the skeleton pixel point from the end point; the nodes are branch points or end points; the search path is composed of all skeleton pixel points passing between the end point and the adjacent node in the search process; the length of the search path is the sum of distances between all adjacent skeleton pixel points on the search path;
calculating a nearest Euclidean distance offshore of the neighboring node; the nearest Euclidean distance of the offshore is the shortest Euclidean distance between the adjacent nodes and the non-water body pixel points;
judging whether a search path meeting the pruning condition exists in the current river skeleton image under the current iteration number m;
if yes, deleting the search path meeting the pruning condition in the current river skeleton image to obtain the river skeleton image after the path deletion, adding 1 to the current iteration number m, taking the river skeleton image after the path deletion as the current river skeleton image, and returning to the determined current iteration number m to determine the end point and the branch point in the current river skeleton image; the pruning condition is determined by the search path length corresponding to the end point, the nearest Euclidean distance of the offshore, a preset proportionality coefficient threshold value and a preset absolute length threshold value; the preset proportionality coefficient threshold is a set ratio of the search path length to the nearest Euclidean distance offshore; the preset absolute length threshold is a set absolute length value of the search path;
if not, determining the river skeleton image obtained by deleting the path under the last iteration number m-1 as the pruned river skeleton image;
determining a connectivity matrix of the pruned river skeleton image; the total line number and the total column number of the connectivity matrix are equal to the sum of the number of end points and branch points in the pruned river skeleton image; the connectivity matrix represents the extracted topology information of the river network;
the pruning conditions comprise a first pruning constraint condition and a second pruning constraint condition;
the first pruning constraint condition is DR(h,j)＜r·DB(j) (ii) a Wherein D isR(h, j) represents the length of the search path corresponding to the endpoint h, and j represents the adjacent node of the endpoint h; r represents a preset scaling factor threshold; dB(j) Represents the nearest euclidean distance offshore of the neighboring node j;
the second pruning constraint condition is DR(h, j) < d; where d represents a preset absolute length threshold.
2. The method for extracting topology information of a river network according to claim 1, wherein the extracting a river skeleton in the water body image by using a morphological refinement algorithm to obtain a river skeleton image specifically comprises:
under the current iteration number n, judging whether water body pixel points meeting a first constraint condition or a second constraint condition exist in the current water body image by adopting a morphological thinning algorithm; the first constraint condition and the second constraint condition are both determined by pixel values in an eight-pixel neighborhood around a water body pixel point;
if so, deleting the water body pixel points meeting the first constraint condition or the second constraint condition, updating the water body image, adding 1 to the current iteration number n, taking the updated water body image as the current water body image, and returning to the current iteration number n to judge whether the water body pixel points meeting the first constraint condition or the second constraint condition exist in the current water body image by adopting a morphological thinning algorithm;
and if not, taking the water body image updated under the last iteration number n-1 as a river skeleton image.
3. The method for extracting topology information of river network according to claim 2, wherein the first constraint condition comprises a first constraint sub-condition, a second constraint sub-condition and a third constraint sub-condition; the second constraint includes a first constraint sub-condition, a second constraint sub-condition, and a fourth constraint sub-condition;
the first constraint sub-condition is XH(p) ═ 1; wherein the content of the first and second substances,
Figure FDA0002961356120000031
p represents a water body pixel point; biA neighborhood pixel operator representing pixel p;
Figure FDA0002961356120000032
wherein x is2i-1、x2iAnd x2i+1The method comprises the steps of representing neighborhood pixel values in eight-pixel neighborhood around a pixel p, wherein i is 1,2,3 and 4, 8 pixel points in the neighborhood of a water body pixel point p are arranged in sequence anticlockwise by taking an adjacent water body pixel point on the right side of p as a starting point, and the pixel value of the adjacent pixel point on the right side of p is x1The pixel values of the adjacent water body pixel points which are sequentially arranged anticlockwise are x respectively2，x3......x8，x9Wherein x is9＝x1(ii) a The pixel value of the water body pixel point is represented by a numerical value 1, and the pixel value of the non-water body pixel point is represented by a numerical value 0;
the second constraint sub-condition is that min { n is more than or equal to 21(p),n2(p) } is less than or equal to 3; wherein the content of the first and second substances,
Figure FDA0002961356120000033
v represents a logical OR operation, and when k is 4, x2k+1＝x1；
The third constraint sub-condition is
Figure FDA0002961356120000034
Λ represents a logical and operation,
Figure FDA0002961356120000035
representing a pixel value x8Is negated；
The fourth constraint sub-condition is
Figure FDA0002961356120000036
Figure FDA0002961356120000037
Representing a pixel value x4The inverted value of (d).
4. The method for extracting topology information of a river network according to claim 1, wherein the determining a connectivity matrix of the pruned river skeleton image specifically includes:
determining an initial connectivity matrix of the pruned river skeleton image; the initial connectivity matrix is a w-row and w-column all-zero matrix, and w is the total number of nodes in the pruned river skeleton image; the nodes are end points or branch points;
for each node e in the pruned river skeleton image, searching a node f adjacent to the node e in any direction in the neighborhood of the node e by adopting a region growing method; e and f are the node numbers in the pruned river skeleton image; e is more than or equal to 1 and less than or equal to w, f is more than or equal to 1 and less than or equal to w, and e is not equal to f;
judging whether a path connecting the node e and the node f meets a path connection condition or not; the path connection condition is that all the pixels on the path are skeleton pixels, and nodes except the node e and the node f do not exist on the path;
and if so, setting the element of the e-th row and the f-th column in the initial connectivity matrix as 1.
5. A system for extracting topology information of a river network, comprising:
the water body image extraction module is used for extracting the water body in the remote sensing image to obtain a water body image; the water body image consists of a plurality of water body pixel points and non-water body pixel points;
the skeleton image extraction module is used for extracting a river skeleton in the water body image by adopting a morphological thinning algorithm to obtain a river skeleton image; the river skeleton image is composed of a plurality of skeleton pixel points;
the first determining module is used for determining end points and branch points in the current river skeleton image under the current iteration number m; only one pixel point in the neighborhood of eight pixels around the end point is a skeleton pixel point; more than two pixel points in the neighborhood of eight pixels around the branch point are skeleton pixel points;
the path searching module is used for respectively taking each end point as a starting point and searching adjacent nodes of each end point along a skeleton path by adopting a path searching algorithm to obtain a searching path and a searching path length corresponding to the end point; the adjacent node is a first node which is searched along the skeleton pixel point from the end point; the nodes are branch points or end points; the search path is composed of all skeleton pixel points passing between the end point and the adjacent node in the search process; the length of the search path is the sum of distances between all adjacent skeleton pixel points on the search path;
an offshore euclidean distance calculation module for calculating the nearest offshore euclidean distance of the neighboring nodes; the nearest Euclidean distance of the offshore is the shortest Euclidean distance between the adjacent nodes and the non-water body pixel points;
the judging module is used for judging whether a search path meeting the pruning condition exists in the current river skeleton image under the current iteration number m;
the pruning module is used for deleting the search path meeting the pruning condition in the current river skeleton image if the search path meeting the pruning condition exists in the current river skeleton image under the previous iteration number m to obtain the river skeleton image after the path deletion, adding 1 to the current iteration number m, taking the river skeleton image after the path deletion as the current river skeleton image, and returning to the first determining module; the pruning condition is determined by the search path length corresponding to the end point, the nearest Euclidean distance of the offshore, a preset proportionality coefficient threshold value and a preset absolute length threshold value; the preset proportionality coefficient threshold is a set ratio of the search path length to the nearest Euclidean distance offshore; the preset absolute length threshold is a set absolute length value of the search path;
the second determining module is used for determining the river skeleton image obtained by deleting the path under the previous iteration number m-1 as the pruned river skeleton image if the current river skeleton image does not have the search path meeting the pruning condition under the previous iteration number m;
the third determination module is used for determining a connectivity matrix of the pruned river skeleton image; the total line number and the total column number of the connectivity matrix are equal to the sum of the number of end points and branch points in the pruned river skeleton image; the connectivity matrix represents the extracted topology information of the river network;
the pruning conditions comprise a first pruning constraint condition and a second pruning constraint condition;
the first pruning constraint condition is DR(h,j)＜r·DB(j) (ii) a Wherein D isR(h, j) represents the length of the search path corresponding to the endpoint h, and j represents the adjacent node of the endpoint h; r represents a preset scaling factor threshold; dB(j) Represents the nearest euclidean distance offshore of the neighboring node j;
the second pruning constraint condition is DR(h, j) < d; where d represents a preset absolute length threshold.
6. The river network topology information extraction system according to claim 5, wherein the skeleton image extraction module specifically includes:
the first judgment unit is used for judging whether water body pixel points meeting a first constraint condition or a second constraint condition exist in the current water body image by adopting a morphological thinning algorithm under the current iteration number n; the first constraint condition and the second constraint condition are both determined by pixel values in an eight-pixel neighborhood around a water body pixel point;
the water body image updating unit is used for deleting the water body pixel points meeting the first constraint condition or the second constraint condition if the water body pixel points meeting the first constraint condition or the second constraint condition exist in the current water body image, updating the water body image, adding 1 to the current iteration number n, taking the updated water body image as the current water body image, and returning the current water body image to the first judging unit;
and the skeleton image determining unit is used for taking the water body image updated under the last iteration number n-1 as a river skeleton image if no water body pixel points meeting the first constraint condition or the second constraint condition exist in the current water body image.
7. The river network topology information extraction system according to claim 6, wherein the first constraint condition includes a first constraint sub-condition, a second constraint sub-condition and a third constraint sub-condition; the second constraint includes a first constraint sub-condition, a second constraint sub-condition, and a fourth constraint sub-condition;
the first constraint sub-condition is XH(p) ═ 1; wherein the content of the first and second substances,
Figure FDA0002961356120000061
p represents a water body pixel point; biA neighborhood pixel operator representing pixel p;
Figure FDA0002961356120000062
wherein x is2i-1、x2iAnd x2i+1The method comprises the steps of representing neighborhood pixel values in eight-pixel neighborhood around a pixel p, wherein i is 1,2,3 and 4, 8 pixel points in the neighborhood of a water body pixel point p are arranged in sequence anticlockwise by taking an adjacent water body pixel point on the right side of p as a starting point, and the pixel value of the adjacent pixel point on the right side of p is x1The pixel values of the adjacent water body pixel points which are sequentially arranged anticlockwise are x respectively2，x3......x8，x9Wherein x is9＝x1(ii) a The pixel value of the water body pixel point is represented by a numerical value 1, and the pixel value of the non-water body pixel point is represented by a numerical value 0;
the second constraint sub-condition is that min { n is more than or equal to 21(p),n2(p) } is less than or equal to 3; wherein the content of the first and second substances,
Figure FDA0002961356120000063
v represents a logical OR operation, and when k is 4, x2k+1＝x1；
The third constraint sub-condition is
Figure FDA0002961356120000071
Λ represents a logical and operation,
Figure FDA0002961356120000072
representing a pixel value x8The negation value of (a);
the fourth constraint sub-condition is
Figure FDA0002961356120000073
Figure FDA0002961356120000074
Representing a pixel value x4The inverted value of (d).
8. The river network topology information extraction system according to claim 6, wherein the third determining module specifically includes:
the initial matrix determining unit is used for determining an initial connectivity matrix of the pruned river skeleton image; the initial connectivity matrix is a w-row and w-column all-zero matrix, and w is the total number of nodes in the pruned river skeleton image; the nodes are end points or branch points;
the searching unit is used for searching a node f adjacent to the node e in any direction in the neighborhood of the node e by adopting a region growing method for each node e in the pruned river skeleton image; e and f are the node numbers in the pruned river skeleton image; e is more than or equal to 1 and less than or equal to w, f is more than or equal to 1 and less than or equal to w, and e is not equal to f;
the second judging unit is used for judging whether a path connecting the node e and the node f meets the path connecting condition or not; the path connection condition is that all the pixels on the path are skeleton pixels, and nodes except the node e and the node f do not exist on the path;
and the communication matrix determining unit is used for setting the element of the e-th row and the f-th column in the initial communication matrix as 1 if the path connecting the node e and the node f meets the path connecting condition.
CN201910937320.2A 2019-09-30 2019-09-30 Method and system for extracting topology information of river network Active CN110688961B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910937320.2A CN110688961B (en) 2019-09-30 2019-09-30 Method and system for extracting topology information of river network

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910937320.2A CN110688961B (en) 2019-09-30 2019-09-30 Method and system for extracting topology information of river network

Publications (2)

Publication Number Publication Date
CN110688961A CN110688961A (en) 2020-01-14
CN110688961B true CN110688961B (en) 2021-06-25

Family

ID=69111235

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910937320.2A Active CN110688961B (en) 2019-09-30 2019-09-30 Method and system for extracting topology information of river network

Country Status (1)

Country Link
CN (1) CN110688961B (en)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112258570B (en) * 2020-11-17 2022-12-20 中国科学院空天信息创新研究院 Method and system for measuring and calculating water surface width of river full water area
CN113989653B (en) * 2021-09-16 2024-04-09 黄河水利委员会黄河水利科学研究院 Method for extracting river plane geometric index and topological structure
CN114140466B (en) * 2022-02-07 2022-05-31 浙江托普云农科技股份有限公司 Plant root system measuring method, system and device based on image processing
CN117315501B (en) * 2023-10-23 2024-04-12 中国水利水电科学研究院 Remote sensing water body classification method based on water body plaque shape and adjacent relation
CN117788395A (en) * 2023-12-15 2024-03-29 南京林业大学 System and method for extracting root phenotype parameters of pinus massoniana seedlings based on images
CN117437309B (en) * 2023-12-21 2024-03-22 梁山公用水务有限公司 Water conservancy water affair digital management system based on artificial intelligence

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103400151A (en) * 2013-08-16 2013-11-20 武汉大学 Optical remote-sensing image, GIS automatic registration and water body extraction integrated method
CN105719300A (en) * 2016-01-22 2016-06-29 黄河水利委员会信息中心 Riverway main stream line detection method based on SNE manifold learning
CN105740807A (en) * 2016-01-28 2016-07-06 武汉大学 Method for extracting continuous river skeleton lines from remote sensing image based on mathematical morphology
CN106295890A (en) * 2016-08-12 2017-01-04 中国水利水电科学研究院 The connective method of river network is evaluated based on Complex Networks Theory
CN107657618A (en) * 2017-10-10 2018-02-02 中国科学院南京地理与湖泊研究所 Regional scale erosion groove extraction method based on remote sensing image and terrain data
WO2018042208A1 (en) * 2016-09-05 2018-03-08 Xihelm Limited Street asset mapping

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101976347A (en) * 2010-10-21 2011-02-16 西北工业大学 Method for recognizing overwater bridge in remote sensing image on basis of Mean Shift segmentation

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103400151A (en) * 2013-08-16 2013-11-20 武汉大学 Optical remote-sensing image, GIS automatic registration and water body extraction integrated method
CN105719300A (en) * 2016-01-22 2016-06-29 黄河水利委员会信息中心 Riverway main stream line detection method based on SNE manifold learning
CN105740807A (en) * 2016-01-28 2016-07-06 武汉大学 Method for extracting continuous river skeleton lines from remote sensing image based on mathematical morphology
CN106295890A (en) * 2016-08-12 2017-01-04 中国水利水电科学研究院 The connective method of river network is evaluated based on Complex Networks Theory
WO2018042208A1 (en) * 2016-09-05 2018-03-08 Xihelm Limited Street asset mapping
CN107657618A (en) * 2017-10-10 2018-02-02 中国科学院南京地理与湖泊研究所 Regional scale erosion groove extraction method based on remote sensing image and terrain data

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
"基于主曲线的遥感图像河岸线提取";郭芸等;《通信学报》;20161130;第37卷(第11期);第80-89页 *
"数学形态学与拓扑约束支持的单条河流骨架线提取";孟令奎等;《遥感学报》;20170531;第21卷(第05期);第785-795页 *

Also Published As

Publication number Publication date
CN110688961A (en) 2020-01-14

Similar Documents

Publication Publication Date Title
CN110688961B (en) Method and system for extracting topology information of river network
CN106981092B (en) Priority-Flood-based internal flow domain extraction method
CN110008602B (en) Road network selection method considering multi-feature coordination under large scale
CN110309248B (en) Method for automatically dividing traffic cells of traffic road network based on Voronoi diagram
CN113297429B (en) Social network link prediction method based on neural network architecture search
CN113822832A (en) Natural resource multi-source vector data fusion method
CN116704137B (en) Reverse modeling method for point cloud deep learning of offshore oil drilling platform
CN111612750A (en) Overlapping chromosome segmentation network based on multi-scale feature extraction
CN112418049A (en) Water body change detection method based on high-resolution remote sensing image
CN111143588A (en) Image space-time index quick retrieval method based on machine learning
CN111144487B (en) Method for establishing and updating remote sensing image sample library
CN114140466B (en) Plant root system measuring method, system and device based on image processing
CN116415318B (en) Modeling method for hydrologic connectivity of lake in inner flow area based on mathematical morphology
CN116431964A (en) Run-length stripping method for generating complex river network water system skeleton line
CN113886999B (en) Method for quickly establishing drainage model and GIS (geographic information System) data through CAD (computer aided design) data
CN115795947A (en) River channel reverse tracing determination method based on digital elevation model
CN115688247A (en) Loess plateau channel area extraction and calculation method
CN114596490A (en) Hilly land feature line extraction method and hilly land DEM (digital elevation model) fine production method
CN111553272B (en) High-resolution satellite optical remote sensing image building change detection method based on deep learning
CN114359558A (en) Roof image segmentation method based on hybrid framework
CN110659823B (en) Similar basin analysis method, model, system and computer storage medium
Cheng et al. Automated detection of impervious surfaces using night-time light and Landsat images based on an iterative classification framework
CN111914039B (en) Road network updating method and device
CN113327327A (en) Automatic extraction method of planar water system skeleton line
CN111985366A (en) Road center line and pile number identification method and device

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant