CN116415318A - 一种基于数学形态学的内流区湖泊水文连通性建模方法 - Google Patents

一种基于数学形态学的内流区湖泊水文连通性建模方法 Download PDF

Info

Publication number
CN116415318A
CN116415318A CN202310283430.8A CN202310283430A CN116415318A CN 116415318 A CN116415318 A CN 116415318A CN 202310283430 A CN202310283430 A CN 202310283430A CN 116415318 A CN116415318 A CN 116415318A
Authority
CN
China
Prior art keywords
lake
overflow
pixels
depression
depressions
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.)
Granted
Application number
CN202310283430.8A
Other languages
English (en)
Other versions
CN116415318B (zh
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.)
Huazhong University of Science and Technology
Changjiang Institute of Survey Planning Design and Research Co Ltd
Changjiang River Scientific Research Institute Changjiang Water Resources Commission
Original Assignee
Huazhong University of Science and Technology
Changjiang Institute of Survey Planning Design and Research Co Ltd
Changjiang River Scientific Research Institute Changjiang Water Resources Commission
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 Huazhong University of Science and Technology, Changjiang Institute of Survey Planning Design and Research Co Ltd, Changjiang River Scientific Research Institute Changjiang Water Resources Commission filed Critical Huazhong University of Science and Technology
Priority to CN202310283430.8A priority Critical patent/CN116415318B/zh
Publication of CN116415318A publication Critical patent/CN116415318A/zh
Application granted granted Critical
Publication of CN116415318B publication Critical patent/CN116415318B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • 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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Geometry (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Data Mining & Analysis (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Image Processing (AREA)

Abstract

一种基于数学形态学的内流区湖泊水文连通性建模方法，包括：获取研究区的数字高程模型；对数字高程模型进行虚假洼地填充，获得无虚假洼地数字高程模型；对无虚假洼地数字高程模型进行平地检测，获得数字高程模型平地像素集，计算各湖泊的湖心点像元；基于各湖泊的湖心点像元进行数字高程模型分水岭变换，并对数字高程模型分水岭进行分割，计算获得湖泊影响区域骨架线以获得潜在溢出口位置；根据潜在溢出口位置构建漫溢邻接矩阵；建立湖泊漫溢级联结构，多个湖泊漫溢级联结构形成漫溢级联森林结构。本发明较好地解决了常规Priority‑flood算法在内流区往往容易出现错误处理水流方向的难题，有利于内流区湖泊的水文连通性建模与分析。

Description

一种基于数学形态学的内流区湖泊水文连通性建模方法
技术领域
本发明涉及遥感及水文领域，具体是一种基于数学形态学的内流区湖泊水文连通性建模方法。
背景技术
湖泊一般由相对封闭可蓄水的洼地形成，洼地是地球上普遍存在的地貌单元，大量地理景观由不同尺度的洼地构成。洼地因地理位置、表面积、深度、蓄水量、几何形状和其他特性的差异而呈现不同的水文和生态效应。量化和理解洼地及其结构可以促进对湿地保护、田间持水量和洪水淹没制图的理解。同时，地表洼地的识别、描绘和量化可以为许多环境应用提供关键信息，例如湿地保护和管理、地表水制图、水文连通性分析和流域建模等。
栅格DEM中的洼地是指没有较低邻域的局部高程极小值的格网。DEM中的洼地分为虚假洼地和真实的自然洼地。早期的栅格DEM因噪声、误差及空间分辨率不足存在了大量虚假洼地，在水文建模和数字地形分析中，这些虚假洼地被当作缺陷而需要进行填充从而创建无洼地的DEM。消除洼地的方法有填充、刻录(Carve)以及组合的方法，这也是DEM预处理的主要功能之一。这类方法会破坏自然洼地的复杂地形结构，改变了真实地形，不利于水文连通性分析和水文建模应用。
洼地处理在水文建模过程中处于非常重要的环节。到目前为止，对洼地处理的已经存在很多的研究，如利用某些阈值(例如，大小，深度，体积)将实际的洼地与虚假的洼地分开。然而，人们很少关注地表洼地的几何和拓扑性质，也很少关注地表洼地对局部和区域水文生态的影响。少有研究试图显式地描绘和量化实际洼地的漫溢级联结构，这可以为表征地表水文连通性和模拟充水-溢流水文过程提供关键信息。因此，研究洼地漫溢级联结构，对了解整个区域水文特征具有非常重要的意义，对认识区域湖泊漫溢过程具有很好的指导作用，对评估区域内湖泊的漫溢风险，减少因湖泊漫溢给人类生命财产造成的危害具有非常重要的意义。
发明内容
针对目前人们很少关注地表洼地的几何和拓扑性质以及地表洼地对局部和区域水文生态的影响，导致显式地描绘和量化实际洼地的漫溢级联结构研究不足的问题，本发明提出一种基于数学形态学的内流区湖泊水文连通性建模方法。在该方法中，创新采用测地数学形态学方法改进Priority-flood算法，从DEM数据中识别提取内流区湖泊及其分水线，采用洼地级联关系模型构建湖泊漫溢级联拓扑关系，以森林结构图的方式揭示内流区湖泊群之间的水文连通性，较好地解决了常规Priority-flood算法在内流区往往容易出现错误处理水流方向的难题，有利于内流区湖泊的水文连通性建模与分析。
为实现上述目的，本发明提出一种基于数学形态学的内流区湖泊水文连通性建模方法，包括以下步骤：
步骤S1：获取研究区的数字高程模型；
步骤S2：利用数学形态学方法中的测地重建方法对步骤S1获取的数字高程模型进行虚假洼地填充，获得无虚假洼地数字高程模型；
步骤S3：利用数学形态学方法中的测地区域生长方法对所述无虚假洼地数字高程模型进行平地检测，获得数字高程模型平地像素集，即湖泊集，通过计算各湖泊的几何中心，得到各湖泊的湖心点像元；
步骤S4：基于步骤S3所得各湖泊的湖心点像元，利用数学形态学方法中淹没模拟方法进行数字高程模型分水岭变换，并利用标记控制方法对数字高程模型分水岭进行分割，计算获得湖泊影响区域骨架线，从而获得潜在溢出口位置；
步骤S5：根据所述潜在溢出口位置构建漫溢邻接矩阵；
步骤S6：基于改进的Priority-flood算法以及步骤S5所得漫溢邻接矩阵，建立湖泊漫溢级联结构，多个湖泊漫溢级联结构形成漫溢级联森林结构。
进一步的，步骤S2中的测地重建方法是指有界图像经过有限次数的测地变换迭代计算后总会收敛，直到标记图像的扩张或者收缩完全被掩膜图像阻止，如果继续循环迭代，标记图像的像素值不再发生改变，测地重建方法包括膨胀重建和腐蚀重建。
进一步的，所述腐蚀重建用于消除洼地，具体过程为：从标记DEM非边界像元开始，按从左到右、从上到下的顺序遍历每一个格网，依次进行测地腐蚀计算，将标记DEM中心格网的8邻域的最小值与原始DEM对应网格高程进行比较，取两者的较大值赋值给中心网格，当全部格网遍历完成后则测地腐蚀重建结束。
进一步的，所述步骤S3中利用数学形态学方法中的测地区域生长方法进行数字高程模型平地检测，获得数字高程模型平地像素集，具体包括：
基于区域生长的平地检测算法引入两个栈数据结构，分别记为A和B，其中栈A存储种子像元，栈B存储目前已确定的连通区域像元，以及一个数组C存储种子像元的邻域像元，采用一个与DEM同维度的二维标记矩阵，记录每个像元的处理状态，基于区域生长的平地检测算法遵循从左到右、从上到下的DEM像元遍历的原则，并且从DEM的左上角像元开始计算；
具体步骤为：
步骤3.1：选择DEM中尚未处理的像元，将其同时压入栈A、栈B进行初始化，在标记矩阵中设置该像元为已处理状态；
步骤3.2：若栈A不为空，从栈A弹出一个像元作为种子点；否则，转向步骤3.4；
步骤3.3：采用广度优先搜索算法，以种子点为中心，按顺时针方向逐一遍历该种子点8邻域内的像元，并将尚未处理的邻域像元存入数组C中；
步骤3.4：若栈A为空，判断栈B中的像元数量是否小于给定阈值，若满足，表示没有找到平地；否则表示找到了平地，则标记栈B中的所有像元为平地，并给这些像元赋予唯一的平地编号；
步骤3.5：清空栈B，转向步骤3.1。
进一步的，步骤3.3具体包括：
步骤3.3-1：判断C是否为空；
步骤3.3-2：如果C不为空时，弹出C中一个邻域像元，如果该邻域像元与种子像元的高程值相等，则将该像元同时压入栈A和栈B，并在标记矩阵中设置这个邻域像元为已处理状态；否则，转向步骤3.3-1；
步骤3.3-3：如果C为空，则表示种子点的8邻域遍历结束，转向执行步骤3.2。
进一步的，所述步骤S4基于步骤S3所得各湖泊的湖心点像元，利用数学形态学方法中淹没模拟方法进行数字高程模型分水岭变换，并利用标记控制方法对数字高程模型分水岭进行分割，计算获得湖泊影响区域骨架线，从而获得潜在溢出口位置，具体包括：
步骤4.1：对步骤S3中获得的湖泊的湖心点像元用唯一标识符进行标识，并将这些标记的湖心像元及DEM边界像元一并存入优先队列中；
步骤4.2：利用Priority-flood分水岭变换方法，对步骤4.1中优先队列中的每个像元按照溢出高程最小的优先原则进行像元的扩展蔓延，逐一确定每个DEM像元所属的湖泊影响区域；
步骤4.3：分水岭变换后，对DEM像元进行遍历，若某像元的8邻域中检测到高程低于该像元且具有不同湖泊标识符号的像元时，则该像元划分为分水岭像元，将分水岭像元连接起来后即获得湖泊影响区域骨架线；采用二维矩阵来存储溢出点，通过遍历集水区相邻湖泊之间的分水岭像元，当二维矩阵尚未记录该分水岭像元时，将该分水岭像元保存为潜在溢出点。
进一步的，所述步骤S5根据潜在溢出口位置，构建漫溢邻接矩阵，具体包括：
步骤5.1：将整个DEM区域划分的各个分水岭赋予唯一标识，采用二维矩阵存储分水岭邻接关系，格网数值记录其潜在溢出高程，在此基础上构建漫溢邻接矩阵；
步骤5.2：设定漫溢邻接矩阵维度为：(n+1)×(n+1)，其中n是划分的集水区数量，其行数和列数相等，且都包含“O”，将该矩阵的第1行和第1列都设置为“O”，格网坐标(I，J)表示分水岭I，J的潜在溢出高程，即有：Value(I，J)＝[I，J]，同一洼地与自身以及并不相邻的洼地之间都没有漫溢关系，其潜在溢出高程值用空值NODATA表示；
漫溢邻接矩阵的存储方式与格网DEM一致。
进一步的，所述步骤S6基于改进的Priority-flood算法以及步骤S5所得漫溢邻接矩阵，建立湖泊漫溢级联结构，多个湖泊漫溢级联结构形成漫溢级联森林结构，具体包括：
步骤6.1：记漫溢邻接矩阵为M，先将M的上三角非空元素压入优先队列pq，将潜在溢出高程越低的格网赋予更高的优先级，设置与M同维度的标记矩阵flag并初始化所有元素为未处理状符“0”，初始化二叉树bst为空，建立辅助函数用于产生新融合洼地全局标识符ID；
步骤6.2：弹出优先队列pq中顶部格网M[i][j]，判断格网坐标行列号是否为0，同时判断对应标记矩阵中的位置是否被处理，如未处理则在标记矩阵flag修改网格值为1，表明该格网为已处理，由于矩阵是对称阵，从判断i行开始；
步骤6.3：如果i不等于0，开始追踪融合型洼地级联树，令leftID＝i，rightID＝j；步骤6.3具体包括：
步骤6.3-1：利用步骤6.1中的辅助函数产生一个新的洼地ID，称之为parentID，以leftID、rightID对应的洼地分别为左右叶子节点、parentID为父节点更新二叉子树bst，设置ParentID＝leftID；
步骤6.3-2：以leftID为对象，查找其邻接洼地，当满足leftID的潜在溢出高程等于其邻接洼地的潜在溢出高程时，则找到了对应的可融合的洼地，记其为rightID，跳转到步骤6.3-1；当leftID的潜在溢出高程不等于其邻接洼地的潜在溢出高程，查找所有邻接洼地中溢出高程最小的洼地，记为metaID，则追踪到了融合洼地的出水口O，以leftID为子结点，metaID为父节点更新二叉子树bst，连接metaID和O，则这棵融合型洼地级联树构建完成；
步骤6.：如果i＝0，开始追踪瀑布型洼地级联树，步骤6.4具体包括：
步骤6.4-1：以出水口O为子节点、以j对应的洼地为parentID，更新只有一个叶子节点的二叉子树bst(退化为链表)；
步骤6.4-2：以parentID为对象，查找其对应洼地的邻接洼地，当邻接洼地的潜在溢出高程高于parentID的潜在溢出高程时，开始追踪parentID的流入洼地inID，以parentID为叶子节点、流入洼地inID为父节点更新二叉子树bst，然后将流入洼地inID设置为parentID；
步骤6.-3：重复步骤6.4-2，直到没有符合条件的洼地，上游湖泊追踪完毕，则这棵瀑布型洼地级联树构建完成；
步骤6.5：返回步骤6.2重复操作，直到优先队列pq为空或者全部格网都参加了运算，计算结束，输出所有生成的二叉树，即湖泊级联关系，将所有的二叉树采用可视化图形绘制方式便形成湖泊漫溢级联森林结构图。
本发明提出了一种基于数学形态学的内流区湖泊水文连通性建模新方法，该方法采用测地数学形态学方法改进Priority-flood算法，从DEM数据中识别提取内流区湖泊及其分水线，采用洼地级联关系模型构建湖泊漫溢级联拓扑关系，以森林结构图的方式揭示内流区湖泊群之间的水文连通性，较好地解决了常规Priority-flood算法在内流区往往容易错误处理水流方向的难题，为类似内流区湖泊的水文连通性建模与分析，提供了一种简单高效、易于编程实现、可视化程度高的新方法。
附图说明
图1是本发明基于数学形态学的内流区湖泊水文连通性建模方法的流程示意图；
图2是湖泊群漫溢级联关系示意图；
图3是对应于图2中的湖泊群的分水岭、潜在溢出口、溢出高程信息示意图；
图4是湖泊群图2构建的漫溢邻接矩阵示意图；
图5是一维视图下相邻湖泊(洼地)M与N的分水岭及水流走向关系示意图；
图6是基于测地腐蚀重建的DEM填洼过程示意图；
图7是研究区DEM数据、湖泊集和湖心点示意图；
图8是一维信号的形态学过分割示意图；
图9是基于湖泊标记的DEM形态学分割示意图；
图10是一维信号的Priority-flood算法改进前后比较示意图；
图11是湖泊分水岭及潜在溢出点分布示意图；
图12是洼地级联关系构建过程示意图；
图13是青藏高原内流湖漫溢级联结构图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
为方便描述湖泊级联关系，本发明首先做如下定义：
定义1：内洼地与外洼地，内洼地是指不与DEM边界和海洋相连的洼地，外洼地是指位于DEM边界或与海洋直接相连的洼地。洼地流出DEM边界之外后最终汇入海洋，因此用大写字母“O”特指海洋(Ocean)或能通过DEM边界流向海洋的洼地。如图1所示，C，E是内洼地，A，B，D，F，G，H，I是外流洼地。
定义2：洼地的邻居，是指与当前洼地毗邻的洼地，两者以分水线相隔。洼地M的所有邻居洼地集合记为
Figure BDA0004138893970000041
如图1所示，/>
Figure BDA0004138893970000042
定义3：潜在溢出高程，是指相邻两个洼地的分水岭上的最低点，用中括号[]表示。如图3所示，[E，F]＝10，表示洼地E与F的最小溢出高程为10，[G，O]＝11，是指洼地G流出DEM边界的潜在溢出高程是11，而[I，O]＝5是指洼地I汇流进入海洋的高程是5。
定义4：漫溢邻接矩阵(Overflowing Adjacency Matrix，OAM)，记录相邻洼地的潜在溢出高程形成的二维矩阵，如图4所示。
定义5：二元运算符用
Figure BDA0004138893970000043
表示，/>
Figure BDA0004138893970000044
表示两者仅仅相邻而没有直接的水力连接(图5(a))；/>
Figure BDA0004138893970000045
表示两种融合形成一个新生洼地(图5(b))；M→N表示相邻洼地M溢出后流向洼地N(图5(c))。此外，二元运算符/>
Figure BDA0004138893970000046
支持优先级运算符号大括号{.}，中括号[.]，和小括号(.)。
当洼地内水体不断注入使水面抬升，相邻的分水岭是否产生水流流动，取决于两者的结构关系。可分为三种：无交换型、融合型、瀑布型。图5(a)属于无交换型，M与N的水流各自流向外部，不产生水力连接；图5(b)是融合型，当水面超过M与N的分水岭，两者产生了水力连接并产生了新的洼地，水流可以双向流动；图5(c)是瀑布型，表示当水面上升至洼地M的湖盆溢出口后进入洼地N，两者同样产生水力连接，但水流只能单向从M流向N。
Figure BDA0004138893970000051
的判断条件为：M、N潜在溢出高程高于其它分水岭，数学表达为：
Figure BDA0004138893970000052
Figure BDA0004138893970000053
的判断条件为：M、N潜在溢出高程低于其它分水岭，数学表达式为：
Figure BDA0004138893970000054
M→N的判断条件为：M、N潜在溢出高程是与M关联的潜在溢出高程的最小值，但高于与N关联的潜在溢出高程，数学表达式为：
Figure BDA0004138893970000055
如图1所示，本发明实施例提供一种基于数学形态学的内流区湖泊水文连通性建模方法，其特征在于：包括如下步骤：
步骤S1：获取研究区的数字高程模型；
步骤S2：利用数学形态学方法中的测地重建方法进行数字高程模型虚假洼地填充，获得无虚假洼地数字高程模型。
步骤S2中的数学形态学方法是建立在集合论、拓扑学、随机几何、格点理论、非线性偏微分方程的信号处理方法，基本思想是利用携带了方向、大小、色度等先验知识的结构元素(Structural element)对信号进行探测分析，是数字图像处理分析的重要工具。其中，测地重建是指有界图像经过有限次数的测地变换迭代计算后总会收敛，直到标记图像的扩张或者收缩完全被掩膜图像阻止。如果继续循环迭代，标记图像的像素值不再发生改变。测地重建包括膨胀重建和腐蚀重建。
图6描述了DEM腐蚀重建消除洼地的方法，其中图6(a)为原始DEM数据，图6(b)为标记DEM数据，首先查找图6(a)原始DEM数据中高程最大值，然后对图6(a)原始DEM数据进行处理，将除边界外其他网格高程均设置为最大值，图6(b)中灰色区域赋值为65。
图6描述的DEM腐蚀重建消除洼地的方法具体过程为：从标记DEM非边界像元开始，按从左到右、从上到下的顺序遍历每一个格网，依次进行测地腐蚀计算。具体地，将标记DEM中心格网的8邻域的最小值与原始DEM对应网格高程进行比较，取两者的较大值赋值给中心网格，当全部格网遍历完成后则测地腐蚀重建结束。图6(c)至图6(k)展示了整个计算过程。从图上看出，原始DEM中洼地格网C3的高程值由58抬升到59，如图6(k)所示，可见洼地被正确填充形成了无洼地的DEM。
基于测地腐蚀重建的DEM填洼算法的优势在于能够消除虚假洼地，并且只需要遍历一次DEM格网，内存消耗小，算法效率极高。由于遥感测量中对水面处理的问题可能导致真实的洼地地形如湖泊可能会被错误处理，需要对处理结果进行干预以获得正确的填洼结果。此外，洼地填平后造成了新的平坦区域，不利于水流确定算法确定流向，改进方法是在抬升洼地格网时，施加一个极小增量来消除DEM填洼形成的平地。
步骤S3：利用数学形态学方法中测地区域生长方法进行数字高程模型平地检测，获得数字高程模型平地像素集，即湖泊集；测地区域生长(Region growing)方法是形态学中提取满足特定准则的、彼此连通的像素集合起来构成区域的方法，从形态学视角来看，DEM中的平地是指相互连通的具有相同高程值的像元集合区域，区域生长方法利用栈(Stack)来实现，栈是一种先进后出的特殊线性表，只允许在栈顶进行插入和弹出元素，其运算包括压栈(Push)和弹栈(Pop)，压栈是指将DEM中像元压入栈顶，弹栈则是将栈顶元素弹出。
基于区域生长的平地检测算法引入了两个栈数据结构，分别记为A和B，其中栈A存储种子像元，栈B存储目前已确定的连通区域像元，以及一个数组C存储种子像元的邻域像元。除此之外，采用一个与DEM同维度的二维标记矩阵，记录每个像元的处理状态。平地检测算法遵循从左到右、从上到下的DEM像元遍历的原则，并且从DEM的左上角像元开始计算。
具体步骤为：
Step1：选择DEM中尚未处理的像元，将其同时压入栈A、栈B进行初始化，在标记矩阵中设置该像元为已处理状态；
Step2：若栈A不为空，从栈A弹出一个像元作为种子点；否则，转向Step4；
Step3：采用广度优先搜索算法，以种子点为中心，按顺时针方向逐一遍历该种子点8邻域内的像元，并将尚未处理的邻域像元存入数组C中。
Step3-1：判断C是否为空；
Step3-2：如果C不为空时，弹出C中一个邻域像元。如果该邻域像元与种子像元的高程值相等，则将该像元同时压入栈A和栈B，并在标记矩阵中设置这个邻域像元为已处理状态；否则，转向Step3-1；
Step3-3：如果C为空，则表示种子点的8邻域遍历结束，转向执行Step2；
Step4：若栈A为空，判断栈B中的像元数量是否小于给定阈值(默认为1)。若满足前述条件，表示没有找到平地；否则表示找到了平地，则标记栈B中的所有像元为平地，并给这些像元赋予唯一的平地编号。
Step5：清空栈B，转向Step1。
当整个DEM遍历完成后，则基于区域生长的平地检测的运算结束，获得数字高程模型平地像素集，即得到湖泊集。通过计算各湖泊的几何中心，便得到各湖泊的湖心点像元。
图7所示底图为研究区原始DEM数据，通过腐蚀重建消除洼地后，利用区域生长方法对研究区内的平地进行检测，获得相应的平地像素集，即湖泊集，然后计算湖泊集中各湖泊的几何中心，获得各湖泊的湖心点，如图7黑色点所示。
步骤S4：基于步骤S3所得各湖泊的湖心点像元，利用数学形态学方法中淹没模拟方法进行数字高程模型分水岭变换，并利用标记控制方法对数字高程模型分水岭进行分割，计算获得湖泊影响区域骨架线，获得潜在溢出口位置。
传统DEM分水岭变换是基于淹没模拟算法获得，这种方法存在易受到噪声干扰和严重过分割问题，如图8所示。解决过分割的方式有后处理或前处理方式，后处理是在对图像施加分水岭变换之后，依设定准则递归地合并相邻区域，但存在计算耗时长、合并运算终止条件不明确等缺点。前处理通过引入控制符对图像事先标记后再进行变换，合理选择标记能有效地解决过分割问题。
为选取有效的控制标记，本发明实施例利用特征检测从图像中提取对象标记，根据目标对象的属性特点和先验知识，例如图像极值点、平坦区域、同质纹理区域等，赋予不同的控制标记。控制标记划定的区块可以是单一像元或者连通的大片对象区域。将对象标记施加到分割函数中，通过分水岭变换得到目标边界。在图9中，将DEM中的平坦区域作为控制标记，对DEM进行形态学分割得到每个湖泊对应的分水岭分割结果。
如图9所示，通过对平坦区域进行标记，利用淹没模拟方法获得每个湖泊对应的分水岭分割结果。基于淹没模拟的方法进行分水岭变换，是基于3×3的邻域栅格运算确定水流方向，忽视了流域整体地形特征。Priority-flood算法与基于淹没模拟的分水岭变换思想一致，但可从宏观视角提取流域结构特征。Priority-flood算法从DEM边界开始，其实是利用其隐含的流向宏观信息，即高程越低且越在DEM外层的栅格单元成为局部汇点的概率越高，保证了流向在宏观地形的正确性。因此，常规的Priority-flood算法的默认前提是每个DEM内部格网都可以通过非上升的连接路径与DEM边界相通，也即如果DEM内部有水流，则其一定是流向DEM边界。因此，常规的Priority-flood算法适用于外向流域的各种水文分析，而内流区水流并不一定流向DEM边界，如遇到内流区地形特征，采用基于Priority-flood算法进行洼地填充，会将内流区域内的湖泊洼地等真实的地貌被填充，从而产生错误的结果。为解决Priority-flood算法在面对内流区地形时遇到的问题，本发明提出基于标记控制的Priority-flood分水岭分割方法。
所述步骤S4基于步骤S3所得各湖泊的湖心点像元，利用数学形态学方法中淹没模拟方法进行数字高程模型分水岭变换，并利用标记控制方法对数字高程模型分水岭进行分割，计算获得湖泊影响区域骨架线，从而获得潜在溢出口位置，具体包括：
Step1：对步骤S3中获得的湖泊的湖心点像元用唯一标识符进行标识，并将这些标记的湖心像元及DEM边界一并存入队列中；
Step2：利用Priority-flood分水岭变换方法，以湖心种子作为控制标记，按照溢出高程最小的优先原则进行湖心像元的扩展蔓延，逐一确定每个DEM像元所属的湖泊影响区域，；
Step3：分水岭变换后，对DEM像元进行遍历，若某像元的8邻域中检测到高程低于该像元且具有不同湖泊标识符号的像元时，则该像元划分为分水岭像元，将分水岭像元连接起来后即获得湖泊影响区域骨架线；通过遍历集水区相邻湖泊之间的分水岭像元，当集合中尚未记录该分水岭像元时，将该分水岭像元保存为潜在溢出点。尽管可能出现两个毗邻的湖泊之间的分水岭存在多个像元，它们或许具有相等的最小溢出高程值，但只保存第一次遍历的像元作为潜在溢出点。
图10所示为Priority-flood算法改进前后分水岭变换结果。在图10(a)是常规Priority-flood算法，只有边界格网A与M进入优先队列，通过分水岭标记算法，内部洼地被填平，最后被标记为分水岭M，见图10(b)。而在图10(c)是基于改进的方法，它增加了标记控制C和H，因此有四个格网进入优先队列，最终分水岭被标识为3部分，并且内部洼地C和H都得到保留，这样更加符合实际情况，如图10(d)所示。由此可见，基于标记控制的Priority-flood分水岭变换算法能够解决内流区湖泊分水岭提取问题，同时也能应用于外流区湖泊分水岭提取。
图11所示为将DEM边界像元和湖心像元作为输入种子，采用标记控制改进Priority-flood算法提取的DEM分水岭，提取了167个湖泊的集水区，集水区确定后对图11范围内DEM像元进行遍历，当某像元的8邻域中检测到高程低于该像元且具有不同湖泊标识符号的像元时，则该像元划分为分水岭像元，这些像元连接起来后即获得湖泊影响区域骨架线，然后遍历167个集水区相邻湖泊之间的分水岭像元，当集合中尚未记录该像元时，将该像元保存为潜在溢出点，获得482个潜在溢出点，如图11所示。
步骤S5：根据潜在溢出口位置，构建漫溢邻接矩阵。
具体包括：
Step1将整个DEM区域划分的各个分水岭赋予唯一标识，采用二维矩阵存储分水岭邻接关系，格网数值记录其潜在溢出高程，在此基础上构建漫溢邻接矩阵；
Step2：设定漫溢邻接矩阵维度为：(n+1)×(n+1)，其中n是划分的集水区数量，其行数和列数相等，且都包含“O”，将该矩阵的第1行和第1列都设置为“O”，格网坐标(I，J)表示分水岭I，J的潜在溢出高程，即有：Value(I，J)＝[I，J]，同一洼地与自身以及并不相邻的洼地之间都没有漫溢关系，其潜在溢出高程值用空值NODATA表示；
漫溢邻接矩阵的存储方式与格网DEM一致。
根据图11所示482个潜在溢出点，建立了包含行列数为168×168漫溢邻接矩阵。
步骤S6：基于改进的Priority-flood算法以及步骤S5所得漫溢邻接矩阵，建立湖泊漫溢级联结构，提出漫溢级联森林结构。
图4所示为根据图3所示湖泊群的分水岭、潜在溢出口、溢出高程信息构建的湖泊群图2构建的漫溢邻接矩阵，从图中可知漫溢邻接矩阵的存储方式与格网DEM(图6(a)所示)一致，为了构建湖泊漫溢级联森林结构，本发明借鉴Priority-flood算法思想，在漫溢邻接矩阵的基础上，采用优先队列存储所有非NODATA的格网，将潜在溢出高程越低的格网赋予更高的优先级，利用二叉树保存洼地级联模型，并采用可视化图形绘制湖泊漫溢级联关系。
图12洼地级联关系构建过程，其中上半部分为漫溢级联矩阵，下半部分为建立的洼地级联森林结构图，中间灰色格网是优先队列弹出的格网，字母行列中数字符号为产生的新融合洼地的标识符。以图4漫溢邻接矩阵为例，记漫溢邻接矩阵为M，先将M的上三角非空元素压入优先队列pq，设置与M同维度的标记矩阵flag并初始化所有元素为未处理状符“0”，初始化二叉树数组bst为空，建立辅助函数用于产生新融合洼地全局标识符ID。
具体实现过程为：
Step1：记漫溢邻接矩阵为M，先将M的上三角非空元素压入优先队列pq，将潜在溢出高程越低的格网赋予更高的优先级，设置与M同维度的标记矩阵flag并初始化所有元素为未处理状符“0”，初始化二叉树bst为空，建立辅助函数用于产生新融合洼地全局标识符ID。
Step2：弹出优先队列pq中顶部格网M[i][j]，判断格网坐标行列号是否为0，同时判断对应标记矩阵中的位置是否被处理，如未处理则在标记矩阵flag修改网格值为1，表明该格网为已处理，由于矩阵是对称阵，从判断i行开始；
Step3：如果i不等于0，开始追踪融合型洼地级联树。令leftID＝i，rightID＝j；
Step3-1：利用Step1中的辅助函数产生一个新的洼地ID，称之为parentID。以leftID、rightID对应的洼地分别为左右叶子节点、parentID为父节点更新二叉子树bst，设置ParentID＝leftID；
Step3-2：以leftID为对象，查找其邻接洼地，当满足leftID的潜在溢出高程等于其邻接洼地的潜在溢出高程时，则找到了对应的可融合的洼地，记其为rightID，跳转到Step3-1；当leftID的潜在溢出高程不等于其邻接洼地的潜在溢出高程，查找所有邻接洼地中溢出高程最小的洼地，记为metaID，则追踪到了融合洼地的出水口O，以leftID为子结点，metaID为父节点更新二叉子树bst，连接metaID和O，则这棵融合型洼地级联树构建完成。
具体地，当获得优先队列pq中顶部格网对应的i、j后，利用矩阵更新原则对矩阵M进行更新，并将标记矩阵的第j行和第j列标记为已处理。接着，保持i不变，逐步查找矩阵M中j列中最小值，若其同时也是j行中的最小值，则按照前述原则更新M和标记矩阵；否则，迭代结束，融合洼地的计算结束，继续更新M和标记矩阵。
漫溢邻接矩阵更新原则为：在漫溢邻接矩阵M中，用i、j列同一行坐标的最小值格网值填充第i列的相应位置的格网；同理，用i、j行同一列号的最小值格网值填充第i行的相应位置的格网，当遇到空值时直接用另一个数值取代；然后删除第j行、第j列的格网数值。数据更新策略如式(4)。
Figure BDA0004138893970000081
Step4：如果i＝0，开始追踪瀑布型洼地级联树。
Step4-1：以出水口O为子节点、以j对应的洼地为parentID，更新只有一个叶子节点的二叉子树bst(退化为链表)；
Step4-2：以parentID为对象，查找其对应洼地的邻接洼地，当邻接洼地的潜在溢出高程高于parentID的潜在溢出高程时，开始追踪parentID的流入洼地inID，以parentID为叶子节点、流入洼地inID为父节点更新二叉子树bst，然后将流入洼地inID设置为parentID；
Step4-3：重复Step4-2，直到没有符合条件的洼地，上游湖泊追踪完毕，则这棵瀑布型洼地级联树构建完成；
具体地，当pq中弹出元素M[0][j]后，在第j列中查找未被处理的具有最小值的格网，设其格网坐标为(r，j)，其值为M[r][j]，标记(r，j)和(j，r)为已处理状态；然后在第r列中找到未处理的格网最小值M[s][r]，若满足公式(5)，那么s是j的流入洼地，继续查找j的流入洼地，直到没有符合上述条件的洼地存在，或者追踪到的洼地为外部洼地，则这棵树计算完成。
M[r][j]＜M[s][r]and M[s][r]＝min(M[i][s]，i＝0，1，2，...，n} (5)
Step5：返回步骤Step2重复操作，直到优先队列pq为空或者全部格网都参加了运算，计算结束，输出所有生成的二叉树，即湖泊级联关系，将所有的二叉树采用可视化图形绘制方式便形成了湖泊漫溢级联森林结构图。
以图3为例，在图12(a)中，因M[3][3]＝3最小，先弹出格网(C，E)，将C，E两行和两列合并，用对应的最小值逐一格网更新其潜在溢出高程值，接着生成一个全局洼地标识符“1”并修改C为1，删除第5行和第5列(洼地E)的格网数值，然后以C，E为叶子，1为父节点构建二叉树，计算结果见图12(b)。同样地，在第3列中(洼地标识为1)，最小值位于第1行(洼地标识为A)，将洼地A，1融合，并产生新洼地2且修改洼地标识符1为2，更新洼地2的潜在溢出高程并删除洼地A的潜在溢出高程值，以A，1为叶子节点，2为父节点修改二叉树，得到结果见图12(c)。依次地，图12(d)、图12(e)、图12(f)分别增加了洼地B，F，I，洼地融合生成的第1棵二叉树构建完成。
在图12(g)中，最小值[O，H]＝4，则在第8列(洼地H)找到最小值[G，H]＝5，因其满足5＞4，因此有G→→H，见图12(h)；接着在第7列(洼地G)找到最小值[D，G]＝6，因为6＞5，因此有D→G；在图12(i)，第4列(洼地D)中，只有[O，D]＝13且是外部洼地，因此这棵二叉树构建结束。最后剩下未处理的格网为[H，I]＝[I，H]＝7，因为[H，I]＝7＞[O，I]＝5，不满足公式(3)，不能构建二叉树。至此，全部漫溢级联矩阵格网都参与了计算，计算结束，图3的洼地级联构建结果如图12(i)中下半部的森林结构。
图13所示为167个湖泊的漫溢级联结构。每个用数字标注的椭圆结点代表图7中的一个湖泊，实线椭圆结点表示原生的分水岭，对应图11的湖泊及其集水区；虚线结点表示两个湖泊的集水区融合后形成的虚拟分水岭，其标识符大于167。图中“O”指Ocean，表示与海洋连通的外流流域。结点之间通过有向边相连，表示两个集水区存在溢出关系，箭头代表水流溢出方向，有向边的数值表示最小溢出高程。通过结点、边绘制了研究区167个大型湖泊的漫溢级联森林结构图。
从图13中可以看出，研究区湖泊网络子系统包含13棵子树，它们通过1、2、4、13、22、85、110、142、158、201、218、219、220结点与外流流域相通。每棵树包含的结点数差别极大，其中，子树1、85、110、142、158没有子节点且直接与O相连，说明这5棵树都只包含一个内流湖泊。一旦这些湖泊来水量足够时，则位于内外流毗邻区的内流湖将漫溢进入外流区。子树2、4、201、219包含的子节点均超过20个，级联关系错综复杂；相比较而言，子树13、22、218、220的结点均不超过10个，漫溢级联关系相对简单。
从级联结构来看，整个内流区以融合型结构为主，占比超过85％，例如
Figure BDA0004138893970000091
Figure BDA0004138893970000092
瀑布型结构较少，在内流区外缘和内部都有分布，分布在外缘的有1→0，23→22→0，85→0，110→0，142→0，158→O等，分布在内流区内部的有167→146→138，119→131→89等。此外还存在更复杂的结构，一种是多连接的瀑布型结构，即两个或更多的湖泊以瀑布型连接进入同一个湖泊，如子树3→→2和112→→2。另一种称之为复合型结构，特指树中同时包含“融合型+瀑布型”，如子树/>
Figure BDA0004138893970000101
Figure BDA0004138893970000102
这两种复杂的级联结构在研究区内占比较少。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何属于本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (8)

1.一种基于数学形态学的内流区湖泊水文连通性建模方法，其特征在于，包括如下步骤：
步骤S1：获取研究区的数字高程模型；
步骤S2：利用数学形态学方法中的测地重建方法对步骤S1获取的数字高程模型进行虚假洼地填充，获得无虚假洼地数字高程模型；
步骤S3：利用数学形态学方法中的测地区域生长方法对所述无虚假洼地数字高程模型进行平地检测，获得数字高程模型平地像素集，即湖泊集，通过计算各湖泊的几何中心，得到各湖泊的湖心点像元；
步骤S4：基于步骤S3所得各湖泊的湖心点像元，利用数学形态学方法中淹没模拟方法进行数字高程模型分水岭变换，并利用标记控制方法对数字高程模型分水岭进行分割，计算获得湖泊影响区域骨架线，从而获得潜在溢出口位置；
步骤S5：根据所述潜在溢出口位置构建漫溢邻接矩阵；
步骤S6：基于改进的Priority-flood算法以及步骤S5所得漫溢邻接矩阵，建立湖泊漫溢级联结构，多个湖泊漫溢级联结构形成漫溢级联森林结构。
2.如权利要求1所述的基于数学形态学的内流区湖泊水文连通性建模方法，其特征在于：步骤S2中的测地重建方法是指有界图像经过有限次数的测地变换迭代计算后总会收敛，直到标记图像的扩张或者收缩完全被掩膜图像阻止，如果继续循环迭代，标记图像的像素值不再发生改变，测地重建方法包括膨胀重建和腐蚀重建。
3.如权利要求2所述的基于数学形态学的内流区湖泊水文连通性建模方法，其特征在于：所述腐蚀重建用于消除洼地，具体过程为：从标记DEM非边界像元开始，按从左到右、从上到下的顺序遍历每一个格网，依次进行测地腐蚀计算，将标记DEM中心格网的8邻域的最小值与原始DEM对应网格高程进行比较，取两者的较大值赋值给中心网格，当全部格网遍历完成后则测地腐蚀重建结束。
4.如权利要求1所述的基于数学形态学的内流区湖泊水文连通性建模方法，其特征在于：所述步骤S3中利用数学形态学方法中的测地区域生长方法进行数字高程模型平地检测，获得数字高程模型平地像素集，具体包括：
基于区域生长的平地检测算法引入两个栈数据结构，分别记为A和B，其中栈A存储种子像元，栈B存储目前已确定的连通区域像元，以及一个数组C存储种子像元的邻域像元，采用一个与DEM同维度的二维标记矩阵，记录每个像元的处理状态，基于区域生长的平地检测算法遵循从左到右、从上到下的DEM像元遍历的原则，并且从DEM的左上角像元开始计算；
具体步骤为：
步骤3.1：选择DEM中尚未处理的像元，将其同时压入栈A、栈B进行初始化，在标记矩阵中设置该像元为已处理状态；
步骤3.2：若栈A不为空，从栈A弹出一个像元作为种子点；否则，转向步骤3.4；
步骤3.3：采用广度优先搜索算法，以种子点为中心，按顺时针方向逐一遍历该种子点8邻域内的像元，并将尚未处理的邻域像元存入数组C中；
步骤3.4：若栈A为空，判断栈B中的像元数量是否小于给定阈值，若满足，表示没有找到平地；否则表示找到了平地，则标记栈B中的所有像元为平地，并给这些像元赋予唯一的平地编号；
步骤3.5：清空栈B，转向步骤3.1。
5.如权利要求4所述的基于数学形态学的内流区湖泊水文连通性建模方法，其特征在于：步骤3.3具体包括：
步骤3.3-1：判断C是否为空；
步骤3.3-2：如果C不为空时，弹出C中一个邻域像元，如果该邻域像元与种子像元的高程值相等，则将该像元同时压入栈A和栈B，并在标记矩阵中设置这个邻域像元为已处理状态；否则，转向步骤3.3-1；
步骤3.3-3：如果C为空，则表示种子点的8邻域遍历结束，转向执行步骤3.2。
6.如权利要求1所述的基于数学形态学的内流区湖泊水文连通性建模方法，其特征在于：所述步骤S4基于步骤S3所得各湖泊的湖心点像元，利用数学形态学方法中淹没模拟方法进行数字高程模型分水岭变换，并利用标记控制方法对数字高程模型分水岭进行分割，计算获得湖泊影响区域骨架线，从而获得潜在溢出口位置，具体包括：
步骤4.1：对步骤S3中获得的湖泊的湖心点像元用唯一标识符进行标识，并将这些标记的湖心像元及DEM边界像元一并存入优先队列中；
步骤4.2：利用Priority-flood分水岭变换方法，对步骤4.1中优先队列中的每个像元按照溢出高程最小的优先原则进行像元的扩展蔓延，逐一确定每个DEM像元所属的湖泊影响区域；
步骤4.3：分水岭变换后，对DEM像元进行遍历，若某像元的8邻域中检测到高程低于该像元且具有不同湖泊标识符号的像元时，则该像元划分为分水岭像元，将分水岭像元连接起来后即获得湖泊影响区域骨架线；采用二维矩阵来存储溢出点，通过遍历集水区相邻湖泊之间的分水岭像元，当二维矩阵尚未记录该分水岭像元时，将该分水岭像元保存为潜在溢出点。
7.如权利要求6所述的基于数学形态学的内流区湖泊水文连通性建模方法，其特征在于：所述步骤S5根据潜在溢出口位置，构建漫溢邻接矩阵，具体包括：
步骤5.1：将整个DEM区域划分的各个分水岭赋予唯一标识，采用二维矩阵存储分水岭邻接关系，格网数值记录其潜在溢出高程，在此基础上构建漫溢邻接矩阵；
步骤5.2：设定漫溢邻接矩阵维度为：(n+1)×(n+1)，其中n是划分的集水区数量，其行数和列数相等，且都包含“O”，将该矩阵的第1行和第1列都设置为“O”，格网坐标(I,J)表示分水岭I,J的潜在溢出高程，即有：Value(I,J)＝[I,J]，同一洼地与自身以及并不相邻的洼地之间都没有漫溢关系，其潜在溢出高程值用空值NODATA表示；
漫溢邻接矩阵的存储方式与格网DEM一致。
8.如权利要求7所述的基于数学形态学的内流区湖泊水文连通性建模方法，其特征在于：所述步骤S6基于改进的Priority-flood算法以及步骤S5所得漫溢邻接矩阵，建立湖泊漫溢级联结构，多个湖泊漫溢级联结构形成漫溢级联森林结构，具体包括：
步骤6.1：记漫溢邻接矩阵为M，先将M的上三角非空元素压入优先队列pq，将潜在溢出高程越低的格网赋予更高的优先级，设置与M同维度的标记矩阵flag并初始化所有元素为未处理状符“0”，初始化二叉树bst为空，建立辅助函数用于产生新融合洼地全局标识符ID；
步骤6.2：弹出优先队列pq中顶部格网M[i][j]，判断格网坐标行列号是否为0，同时判断对应标记矩阵中的位置是否被处理，如未处理则在标记矩阵flag修改网格值为1，表明该格网为已处理，由于矩阵是对称阵，从判断i行开始；
步骤6.3：如果i不等于0，开始追踪融合型洼地级联树，令leftID＝i,rightID＝j；步骤6.3具体包括：
步骤6.3-1：利用步骤6.1中的辅助函数产生一个新的洼地ID，称之为parentID，以leftID、rightID对应的洼地分别为左右叶子节点、parentID为父节点更新二叉子树bst，设置ParentID＝leftID；
步骤6.3-2：以leftID为对象，查找其邻接洼地，当满足leftID的潜在溢出高程等于其邻接洼地的潜在溢出高程时，则找到了对应的可融合的洼地，记其为rightID，跳转到步骤6.3-1；当lefrID的潜在溢出高程不等于其邻接洼地的潜在溢出高程，查找所有邻接洼地中溢出高程最小的洼地，记为metaID，则追踪到了融合洼地的出水口O，以leftID为子结点，metaID为父节点更新二叉子树bst，连接metaID和O，则这棵融合型洼地级联树构建完成；
步骤6.：如果i＝0，开始追踪瀑布型洼地级联树，步骤6.4具体包括：
步骤6.4-1：以出水口O为子节点、以j对应的洼地为parentID，更新只有一个叶子节点的二叉子树bst(退化为链表)；
步骤6.4-2：以parentID为对象，查找其对应洼地的邻接洼地，当邻接洼地的潜在溢出高程高于parentID的潜在溢出高程时，开始追踪parentID的流入洼地inID，以parentID为叶子节点、流入洼地inID为父节点更新二叉子树bst，然后将流入洼地inID设置为parentID；
步骤6.-3：重复步骤6.4-2，直到没有符合条件的洼地，上游湖泊追踪完毕，则这棵瀑布型洼地级联树构建完成；
步骤6.5：返回步骤6.2重复操作，直到优先队列pq为空或者全部格网都参加了运算，计算结束，输出所有生成的二叉树，即湖泊级联关系，将所有的二叉树采用可视化图形绘制方式便形成湖泊漫溢级联森林结构图。
CN202310283430.8A 2023-03-20 2023-03-20 一种基于数学形态学的内流区湖泊水文连通性建模方法 Active CN116415318B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310283430.8A CN116415318B (zh) 2023-03-20 2023-03-20 一种基于数学形态学的内流区湖泊水文连通性建模方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310283430.8A CN116415318B (zh) 2023-03-20 2023-03-20 一种基于数学形态学的内流区湖泊水文连通性建模方法

Publications (2)

Publication Number Publication Date
CN116415318A true CN116415318A (zh) 2023-07-11
CN116415318B CN116415318B (zh) 2024-03-08

Family

ID=87052550

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310283430.8A Active CN116415318B (zh) 2023-03-20 2023-03-20 一种基于数学形态学的内流区湖泊水文连通性建模方法

Country Status (1)

Country Link
CN (1) CN116415318B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117932974A (zh) * 2024-03-21 2024-04-26 威海水利工程集团有限公司 一种水库水下数字高程模型的构建方法

Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010282349A (ja) * 2009-06-03 2010-12-16 Hitachi Engineering & Services Co Ltd 流域抽出プログラムおよび流域抽出装置
CN102842104A (zh) * 2012-07-16 2012-12-26 长江水利委员会长江科学院 面向海量dem数据的高精度河道洪水淹没区生成方法
CN103236086A (zh) * 2013-04-24 2013-08-07 武汉大学 一种顾及地表水文上下文的多尺度dem建模方法
CN105138722A (zh) * 2015-07-14 2015-12-09 南京师范大学 基于数字河湖网络的平原河网区流域集水单元划分方法
CN106981092A (zh) * 2017-03-28 2017-07-25 南京师范大学 基于Priority‑Flood的内流流域提取方法
US20170365094A1 (en) * 2016-04-04 2017-12-21 University Of Cincinnati Localized Contour Tree Method for Deriving Geometric and Topological Properties of Complex Surface Depressions Based on High Resolution Topographical Data
KR101912627B1 (ko) * 2017-05-30 2018-10-30 에스지에이블록체인 주식회사 Gis 기반의 수문-수리모형 분석결과의 통합 가시화 방법
CN110147423A (zh) * 2019-05-21 2019-08-20 中国科学院南京地理与湖泊研究所 基于湖泊汇流关系的内流湖盆区流域自动划分方法
US20190392094A1 (en) * 2018-06-22 2019-12-26 Drexel University Watershed marching-delineation algorithm
AU2020103026A4 (en) * 2020-10-27 2020-12-24 Nanjing Forestry University A Single Tree Crown Segmentation Algorithm Based on Super-pixels and Topological Features in Aerial Images
CN114511509A (zh) * 2022-01-12 2022-05-17 中国石油大学(华东) 基于数字岩心电流场可视化定量表征导电组分贡献的方法
CN115359221A (zh) * 2022-08-17 2022-11-18 中山大学 一种基于数字高程模型的内陆湖泊流域提取方法及系统
CN115470612A (zh) * 2022-07-26 2022-12-13 中国地质大学（北京） 构造湖漫溢溃坝隐患预警方法、系统、介质及设备

Patent Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010282349A (ja) * 2009-06-03 2010-12-16 Hitachi Engineering & Services Co Ltd 流域抽出プログラムおよび流域抽出装置
CN102842104A (zh) * 2012-07-16 2012-12-26 长江水利委员会长江科学院 面向海量dem数据的高精度河道洪水淹没区生成方法
CN103236086A (zh) * 2013-04-24 2013-08-07 武汉大学 一种顾及地表水文上下文的多尺度dem建模方法
CN105138722A (zh) * 2015-07-14 2015-12-09 南京师范大学 基于数字河湖网络的平原河网区流域集水单元划分方法
US20170365094A1 (en) * 2016-04-04 2017-12-21 University Of Cincinnati Localized Contour Tree Method for Deriving Geometric and Topological Properties of Complex Surface Depressions Based on High Resolution Topographical Data
CN106981092A (zh) * 2017-03-28 2017-07-25 南京师范大学 基于Priority‑Flood的内流流域提取方法
KR101912627B1 (ko) * 2017-05-30 2018-10-30 에스지에이블록체인 주식회사 Gis 기반의 수문-수리모형 분석결과의 통합 가시화 방법
US20190392094A1 (en) * 2018-06-22 2019-12-26 Drexel University Watershed marching-delineation algorithm
CN110147423A (zh) * 2019-05-21 2019-08-20 中国科学院南京地理与湖泊研究所 基于湖泊汇流关系的内流湖盆区流域自动划分方法
AU2020103026A4 (en) * 2020-10-27 2020-12-24 Nanjing Forestry University A Single Tree Crown Segmentation Algorithm Based on Super-pixels and Topological Features in Aerial Images
CN114511509A (zh) * 2022-01-12 2022-05-17 中国石油大学(华东) 基于数字岩心电流场可视化定量表征导电组分贡献的方法
CN115470612A (zh) * 2022-07-26 2022-12-13 中国地质大学（北京） 构造湖漫溢溃坝隐患预警方法、系统、介质及设备
CN115359221A (zh) * 2022-08-17 2022-11-18 中山大学 一种基于数字高程模型的内陆湖泊流域提取方法及系统

Non-Patent Citations (8)

* Cited by examiner, † Cited by third party
Title
J. B. LINDSAY等: "Efficient hybrid breaching-filling sink removal methods for flow path enforcement in digital elevation models", HYDROLOGICAL PROCESSES *
R. BARNES等: "Computing water flow through complex landscapes, Part 2: Finding hierarchies in depressions and morphological segmentations", EARTH SURFACE DYNAMICS *
ZHOU G 等: "An efficient variant of the PriorityFlood algorithm for filling depressions in raster digital elevation model", COMPUTERS&GEOSCIENCES, 31 December 2016 (2016-12-31) *
卢庆辉;熊礼阳;蒋如乔;巫晓玲;段家朕;: "一种融合Priority-Flood算法与D8算法特点的河网提取方法", 地理与地理信息科学, no. 04 *
周厚芳;: "流域边界提取方法研究综述", 人民长江, no. 2 *
彭检贵;邢元军;宋亚斌;王宗跃;: "一种顾及地形特征的山区LiDAR高精度DEM提取算法", 软件导刊, no. 10 *
邓必平;严恩萍;洪奕丰;程玉娜;吴群山;: "基于GIS和DEM的东江湖流域水文特征分析", 湖北农业科学, no. 15, 5 August 2013 (2013-08-05) *
陈永刚;杨春菊;陈孝银;孙燕飞;林晨鸷;: "基于树形结构的地形骨架线分级方法研究", 中国矿业大学学报, no. 06 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117932974A (zh) * 2024-03-21 2024-04-26 威海水利工程集团有限公司 一种水库水下数字高程模型的构建方法

Also Published As

Publication number Publication date
CN116415318B (zh) 2024-03-08

Similar Documents

Publication Publication Date Title
WO2024020744A1 (zh) 基于高分辨率影像的生态图斑提取方法
Mackay et al. Extraction and representation of nested catchment areas from digital elevation models in lake‐dominated topography
CN105677890B (zh) 一种城市绿量数字地图制作及显示方法
CN103236086B (zh) 一种顾及地表水文上下文的多尺度dem建模方法
CN106548141B (zh) 一种基于三角网的面向对象耕地信息自动提取方法
CN110415265B (zh) 基于无人机高精度dem数据的梯田自动提取方法
CN105118090B (zh) 一种自适应复杂地形结构的点云滤波方法
CN106981092B (zh) 基于Priority-Flood的内流流域提取方法
Xing et al. A coastal wetlands mapping approach of Yellow River Delta with a hierarchical classification and optimal feature selection framework
CN113807437B (zh) 一种基于dbscan聚类分析的山脊线和山谷线提取方法
Wang et al. The isotropic organization of DEM structure and extraction of valley lines using hexagonal grid
CN113223042A (zh) 一种遥感影像深度学习样本智能采集方法及设备
CN116415318B (zh) 一种基于数学形态学的内流区湖泊水文连通性建模方法
CN109657598B (zh) 基于分层策略的滨海湿地遥感分类方法
Haag et al. Development of a data model to facilitate rapid watershed delineation
CN105279317A (zh) 一种基于dem的平地河网水流方向估算方法
CN105894553A (zh) 一种基于格栅选择的街巷空间形态布局方法
CN112861669A (zh) 地表坡向约束的高分辨率dem地形特征增强提取方法
CN112116709A (zh) 一种提高地形表达精度的地形特征线处理方法
Hao et al. A subpixel mapping method for urban land use by reducing shadow effects
CN115544789B (zh) 基于数字高程模型和坡度成本距离的河谷平原提取方法
CN115359221A (zh) 一种基于数字高程模型的内陆湖泊流域提取方法及系统
CN113128009B (zh) 一种考虑山区平原地貌差异的子流域单元划分方法
CN111274545B (zh) 一种栅格尺度基于地形地貌的多模式产流计算方法
CN110288645B (zh) 一种基于坡向约束的地形表面面积计算方法

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