技术领域本发明属于城市内涝监测预警技术领域，涉及一种基于数字表面模型的城市汇水区提取方法。
背景技术：
数字表面模型(DigitalSurfaceModel，简称DSM)是对地球表面形态的一种离散化的数学表示，也是流域及洪涝等领域进行水文分析的核心数据之一。汇水区，或称集水区域、集水盆地、流域盆地，是指地表径流或其他物质汇聚到一共同出水口的过程中所流经的地表区域。汇水区是水文分析的基本单元，在洪涝分析、数字地形分析、测绘及工程设计中都有着重要的应用。因此从DSM中自动高效地提取汇水区一直都是内涝监测预警技术的一个重要课题。现有的提取方法从采用数据上可以分为两大类：基于栅格和基于不规则三角网(TriangulatedIrregularNetwork，简称TIN)。其中，基于栅格的汇水区提取算法较多，而基于栅格的汇水区提取研究中绝大多数都是基于数字高程模型(DigitalElevationModel，简称DEM)，如基于地表径流漫流模型、基于矩形窗口扫描DEM的模型，以及一些基于水流方向的算法，如多流向算法、DEMON法、Burnin算法、Dinf法、Dco法、Freeman法等，然而这些算法和模型，对城市平原地区有很大的不适用性，Duke等针对城市平原汇水区划分提出了RIDEM模型，左俊杰等在RIDEM进行了改良，但在应用中仍未考虑城区建筑物对汇流的影响，具有一定的局限性，这也是基于DEM提取城市汇水区的难点和局限所在。基于不规则三角网的方法有刘学军提出的基于不规则三角网的河网提取算法，周启鸣提出的基于三角形的DEM模拟了表面水流路径，吴立新提出了一种基于三角形面皮的汇水区划分方法等，基于不规则三角网的方法在描述自然地表形态上有着高精度的优势，但对于城市这种多规则建筑物地区的表达有着一定的缺陷，并不适用。两种数据各有优缺点，栅格数据结构简单，算法处理方便，但采用数字高程模型精度较低，对于一些复杂地形表达效果不够好，使得提取结果不够精细；而不规则三角网能够更加精确的表达真实地形，但却因为数据结构复杂，导致算法处理不便，并没有得到广泛的研究，并且在表达规则建筑物繁多的城市表面形态上存在着一定的缺陷，更难以进行汇水区的提取。综上所述，目前汇水区提取算法中大多基于简单方便的栅格数据，其中基于DEM的占绝大多数，基于数字表面模型数据(DigitalSurfaceModel，简称DSM)的方法极为少见，而城市区域地形复杂，地表建筑繁多，基于DEM的提取方法难以考虑复杂建筑对汇水区的影响，不能满足城区的汇水区域提取要求。而基于TIN的方法在城区表示中存在一定的缺陷，并且太过复杂，研究不多，还没有成熟的解决方案，更不适合大面积城区汇水区域提取。
技术实现要素：
针对现有汇水区提取方法主要针对于地形起伏得到的汇水区域，并不适应于地表建筑繁多，地表与地形相差较大的城市地区汇水区域提取的缺点，本发明的目的是提供一种从城区高精度DSM中提取汇水区域的方法，能够在复杂多变的城区数字表面模型中得到汇水区域。本发明所采用的技术方案是：一种基于数字表面模型的城市汇水区提取方法，其特征在于，包括以下步骤：步骤1：真实流向计算；遍历所有DSM格网点，计算每个格网点的真实流向；步骤2：出水点提取；找出所有的流向指向自身且为局部极小值的点，并进行编号，作为初步的汇水区编号；步骤3：汇水区增长；将出水点作为汇水区域的种子点，将流向指向该汇水区的格网点加入该汇水区点序列，并进行循环迭代，扩大汇水区域范围，直到所有的汇水区面积不再变化；步骤4：汇水区归并；首先对每个汇水区计算最大储水深度和最大储水量，并根据给定的储水深度或者降水量阈值进行汇水区归并；步骤5：汇水区提取；对汇水区边界进行跟踪，得到各个汇水区的矢量进行输出。作为优选，步骤1的具体实现过程包括以下子步骤：步骤1.1：对待计算点e与周围八邻域点的高程进行比较，若周围八邻域所有点的高程均大于或等于e点高程，则e点流向指向自身，否则，e点的流向指向八邻域点中的最低点，若有多个最低点，则将e点的流向暂时设为指向自身；步骤1.2：利用步骤1.1中所述方法对所有的点进行流向判断；步骤1.3：遍历所有流向为(0，0)的格网点，若该点流向指向自身，则判断该点周围八邻域点的流向，然后根据周围八邻域的流向进行矢量合成得到新的流向，直到所有的流向不再变化。作为优选，步骤2中所述提取出水点，若满足以下两个条件则标记为出水点：条件1：该点的流向为(0，0)；条件2：该点为局部极小值的点；对满足条件1，不满足条件2的点，将其流向改为指向其八邻域的任一最低点；对满足上述两个条件的点依次进行编号，作为初步的汇水区编号。作为优选，步骤3的具体实现过程包括以下子步骤：步骤3.1：将所有的出水点作为汇水区增长的种子点，并将出水点的编号作为汇水区的编号，将该出水点加入该汇水区的边界点序列；步骤3.2：按照编号大小依次遍历所有汇水区的边界点序列，依次对该边界点序列的每个点周围的格网点流向进行判断：如果流向指向该边界点，则将该格网点的汇水区编号设置为该边界点汇水区编号，并将该点加入该汇水区的边界点序列，判断完毕后将该边界点从边界点序列删除，然后并将该格网点的汇水区编号改为该出水点的编号；若该点周围点均已有汇水区编号，则将该点从出水点序列删除，然后对下一个汇水区的边界点序列进行同样的操作，直到所有的汇水区边界点序列判断完毕；步骤3.3：对新的边界序列重复上述步骤，直到所有的汇水区边界不再变化，所有汇水区无法增长。作为优选，步骤4的具体实现过程包括以下子步骤：步骤4.1：找到每个汇水区的最低点高程值Cmin和边界序列的最低点高程值Lmin：步骤4.2：求得两者的高程值差作为储水深度D，D＝Cmin-Lmin；步骤4.3：将高程值差小于给定阈值D0的汇水区与该汇水区边界序列最低点Lmin处相邻的另一汇水区合并，并重新计算储水深度D，迭代进行，直到所有汇水区的储水深度D大于给定阈值D0。作为优选，步骤4的具体实现过程包括以下子步骤：步骤4.1：找到每个汇水区的最低点高程值Cmin和边界序列的最低点高程值Lmin：步骤4.2：然后找到汇水区内所有高程值小于或等于Lmin的格网点，并将其与Lmin的差值进行累加，最后将累加值乘以DSM格网点面积即可得到每个汇水区的最大储水量Smax；步骤4.3：通过判断汇水区最大储水量Smax与汇水区域降水量P的关系来决定是否进行汇水区归并；若Smax>P，则不需要归并；若Smax≤P，则将该汇水区与该汇水区边界序列最低点Lmin处相邻的另一汇水区归并；其中：P＝P0*A；其中P为汇水区域降水量，P0为单位面积降水量，A为汇水区面积。与现有的汇水区提取方法相比，本发明的有益效果是：1)本发明与现有方法相比，克服了现有方法针对于地形变化，不适用于地表建筑繁多，地表与地形差异巨大的城市区域的缺点，能够从整体高差起伏剧烈，局部区域平坦的城市区域DSM模型中提取得到汇水区域。2)本发明采用的新的流向计算方法，相比较通常的流向计算方法，考虑的周围格网点的流向对当前格网点流向的影响，并创造性利用矢量合成的方法得到复杂情况流向的判断，适用于城区建筑等复杂表面对水流方向的影响。3)本发明提出的汇水区域提取方法，可以根据积水深度进行汇水区域的合并，能够模拟不同程度降水情况下的汇水区域，满足不同需求下的汇水区提取，具有更强的适应性。附图说明图1为本发明实施例的流程图。图2为本发明实施例的格网点与其8邻域示意图。图3为本发明实施例的流向设置图。具体实施方式为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。请见图1，本发明提供的一种基于数字表面模型的城市汇水区提取方法，包括以下步骤：步骤1：真实流向计算；利用如图2所示的3*3的模板依次对所有DSM格网点进行遍历，首先对中心点e与周围八邻域点(e1～e8)的高程进行比较，按照下面几种情况及流向设置图(图3)进行处理：(1)e1～e8的高程均大于或等于e点高程，则e点流向为(0，0)。(2)e1～e8中存在小于e点高程的点，且只有一个最低点，则e点的流向指向该最低点，具体流向根据图2确定，若e的流向指向e1，则其流向为(1，0)，其他类似确定。(3)e1～e8中存在小于e点高程的点，且有多个高程一样的最低点，则将e点的流向设为(0，0)。待所有的点均按照上面方法进行流向判断之后，则遍历所有流向为(0，0)的格网点，进行流向合成，首先判断该点周围八邻域点的流向，然后根据周围八邻域的流向根据向量的加法法则进行矢量合成，然后计算合成的新的向量与各个流向向量的夹角，取夹角最小的流向向量作为新的流向。重复流向合成步骤，直到所有的流向不再变化。平面向量的加法公式为：设a＝(x1，y1)，b＝(x2，y2)，则：+b＝(x1+x2，y1+y2)；已知平面内两个非零向量，a＝(x1，y1)，b＝(x2，y2)，则两向量之间的夹角α(a，b)计算公式如下：α(a,b)=arccos((x1x2+y1y2)x12+y12*x22+y22);]]>步骤2：出水点提取；采用一个3*3的模板对DSM中每个点进行遍历，若满足以下两个条件则标记为出水点：条件1：该点的流向为(0，0)；条件2：该点为局部极小值的点；对满足条件1，不满足条件2的点，将其流向改为指向其八邻域的任一最低点；对满足上述两个条件的点依次进行编号，作为初步的汇水区编号。步骤3：汇水区增长；将所有的出水点作为汇水区增长的种子点，并将出水点的编号作为汇水区的编号，将该出水点加入该汇水区的边界点序列，按照编号大小依次遍历所有汇水区的边界点序列，依次对该边界点序列的每个点周围的格网点流向进行判断，如果流向指向该边界点，则将该格网点的汇水区编号设置为该边界点汇水区编号，并将该点加入该汇水区的边界点序列，判断完毕后将该边界点从边界点序列删除。然后并将该格网点的的汇水区编号改为该出水点的编号，若该点周围点均已有汇水区编号，则将该点从出水点序列删除，然后对下一个汇水区的边界点序列进行同样的操作，直到所有的汇水区边界点序列判断完毕。然后对新的边界序列重复上述步骤，直到所有的汇水区边界不再变化，所有汇水区无法增长。步骤4：汇水区归并；首先找到每个汇水区的最低点高程值Cmin和边界序列的最低点高程值Lmin，并求得两者的高程值差作为储水深度D，然后将高程值差小于给定阈值D0的汇水区与该汇水区边界序列最低点Lmin处相邻的另一汇水区合并，并重新计算储水深度D，迭代进行，直到所有汇水区的储水深度D大于给定阈值。D＝Cmin-Lmin；若给定单位面积降水量P0，也可通过计算每个汇水区的最大储水量来进行判断是否进行汇水区归并。汇水区的最大储水量的计算方法如下，首先找到每个汇水区的最低点高程值Cmin和边界序列的最低点高程值Lmin，然后找到汇水区内所有高程值小于或等于Lmin的格网点，并将其与Lmin的差值进行累加，最后将累加值乘以DSM格网点面积即可得到每个汇水区的最大储水量Smax，通过判断汇水区最大储水量Smax与汇水区域降水量P的关系来决定是否进行汇水区归并，若Smax>P，则不需要归并，若Smax≤P，则将该汇水区与该汇水区边界序列最低点Lmin处相邻的另一汇水区归并。P＝P0*A；其中P为汇水区域降水量，P0为单位面积降水量，A为汇水区面积。步骤5：汇水区输出。将归并后的汇水区边界进行矢量化并进行输出，得到汇水区矢量。本发明基于高精度的数字表面模型数据，通过真实流向计算及水流模拟分析，提出一种基于数字表面模型的城市汇水区提取方法，将城市建筑物对水流的影响纳入考虑范围，能够从DSM数据中提取得到汇水区域。应当理解的是，本说明书未详细阐述的部分均属于现有技术。应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。