规则高程格网中完整性谷地线的提取.docx
《规则高程格网中完整性谷地线的提取.docx》由会员分享，可在线阅读，更多相关《规则高程格网中完整性谷地线的提取.docx（18页珍藏版）》请在冰豆网上搜索。
规则高程格网中完整性谷地线的提取
规则高程格网中完整性谷地线的提取
王涛
雷蓉
中国测绘科学研究院北京100039
华南农业大学信息学院广州510642
Email：
mapwang@
摘要：
高精度、高分辨率的格网高程数据已经变得非常容易获取，从中自动提取谷地线的工作对水文、地形、生态、土壤等专业领域的科学研究和工程应用有着重要意义。
以往研究分别从局部形态判断、地表流水模拟、地形曲面分析等角度提出大量算法，得到了良好的结果。
然而，现有的主要算法几乎都没有从理论上保证谷地线的完整性，处理过程中使用的图像算法没有顾及地形物理特征。
本文从地形整体特征出发，考虑现有算法在不同地形区域中提取结果的特点，研究设计综合算法模型，扬长避短，发挥不同算法的优点，使提取的谷地线结果具有完整性特点；分析归纳初步结果中的噪声种类和特点，有针对性的设计处理规则，保证系统性谷地线网络的合理性。
关键词：
高程格网谷地线完整性
AbstractValleylinesareabstractrepresentationofterrainandplayimportantrolesinvariousapplication,suchasintheareaofhydrology,biologyandgeomorphometry.Theautomaticextractionofvalleylinehasbeenanactiveresearchtopicforalongtime.Accordingtothedatasourcesused,existingmethodscanbeclassifiedintothreecategories,whicharethosebasedoncontour,basedonTINandbasedongridmatrix.DuetothefastdevelopmentofInSARandLIDARtechnology,highqualitygridelevationdataiswidelyavailableinpresentdays.However,thereisnotamethodwhichcanextractcompletevalleylinesfromthegridelevationmatrix.Andtherearenoteffectivemethodsforreducingtheprocessofnoiseinvalleyline.Thispaperdesignsanovelmodelbasedonanalysisofterraincharacteristicsofvalleylines.Theprimaryvalleylinesareextractedbasedonwateraccumulationalgorithminfirststage.Thenthelocalmorphometryisinvestigatedtoappendeligibleterrainspecificpoint.Thethirdstepconsidersthereductionofnoisesintheresult.Basedontheinductiveanalysis,thenoisesareclassifiedintodifferenttypesandhandledaccordingtothecontext.Discussionandfutureworkaregiveninthelastsection.
Keywordsgridelevationmatrix,valleyline,complete
1.引言
地性线是地表形态的控制元素，是构建结构化地形信息和表达对象化地形要素的关键，也是众多应用中的焦点。
谷地是生物、水文过程的主要发生地，并且它在地形上具有较为明确的物理特征，因此谷地线的自动提取是地性线网络研究中的核心部分，相对而言，谷地线的自动提取研究在成果和发展时间上也最为突出。
当前，已经存在大量的谷地线自动提取算法，它们在算法思想和实现过程上有着很大差异，本文以数据源为标准对这些方法进行了大致分类（如图1）。
图1地性线自动提取方法分类
其中，基于等高线的不同算法在识别特征点的环节上还存在较大差异，通常是使用曲率、三角网、地图代数等方法（费立凡，1993；Thibault等，2000；艾廷华等，2003；吴艳兰，2004），然而这类方法在过滤特征点使用到的阈值需要根据经验设定，在等高线之间连接特征点构建单根谷地线时还存在多义性的问题。
另外当前的数字高程数据源正在发生变化，等高线原有的部分优势将会降低。
从定义和理论上看，三角网高程模型是提取地性线的最佳数据源，因为三角网顶点都应当是由地性点构成，然而由于数据采集方式的限制，当前高程三角网数据基本来自于高程格网和等高线的内插，因此它的应用也受到了很大限制。
随着数据采集技术的发展，规则格网已经成为最主要的高程数据表达形式，同时，规则格网因其数据结构的简单性、计算的并行性，在谷地线的提取研究中成果最为丰富。
众多的算法各自有着独特的优点，然而在理论基础、提取结果等方面分别存在着一定限制。
（1）局部形态法需要对初步结果进行图像细化以降低噪声，而该处理没有顾及谷地线物理特征；其次，基于固定范围检测谷地点导致最终结果限定于固定尺度的地形要素；第三，汇水区信息的获取也使用到没有顾及地形特征的图像处理算法，为后续应用引入了更多的不确定性。
（2）基于流水累积的方法具有合理的物理基础，但是洼地填充和平地流向问题为局部谷地线形态引入了误差，而阈值设定是经验性很强且依赖于局部地貌特征的任务，同时阈值过滤会造成上游谷地线丢失和下游谷地线误判，另外，以往研究都没有处理提取结果中的粘连等噪声。
（3）曲面分析法有着严密的数学基础，能同时标定多种类型的地性点，但是，这种方法在全区域内使用的仍是固定窗口，并不能根据局部地形动态调节窗口大小，同时，检测出的谷地线和山脊线都十分“粗壮”，而且仍然存在零碎线段，因此还需要细化处理并设计合理的连接方案。
现有谷地线提取算法的设计基础都没有保证结果的完整性，另外，使用图像细化算法、平面拓扑连接判断等策略对初步结果进行噪声处理时几乎都没有顾及谷地线的物理特征。
针对这些问题，本文提出，首先利用流水累积法的全局性来提取基本谷地线，然后以第一步结果为启发知识，利用局部形态法识别丢失的典型谷地点，最后对初步结果中的噪声类型进行归纳，分别采取不同策略进行处理。
2.基于流水累积与形态法的谷地线提取
2.1流水累积法
提取谷地线的流水累积法是一种模拟地表径流物理过程的算法，其基本理论基础是，假设每个网格上存在单元量的水，由于地形起伏的存在，地表上的水必然向着低的地方流动并汇集成流，谷地的流水过程远大于坡面过程，这样超过一定汇集量的地方就认为是谷地。
该思想源于1968年Speight以手工的方式从等高线中提取地性线的处理（Mark，1984），Mark模拟这种思想实现了自动化提取谷地线的过程，而后的大量研究促使该思想逐渐完善（O'Callaghan等，1984，Jenson等，1988；Fairfield等，1991；Martz等，1998；Tarboton，1997；Endreny等，2003），成为一类主要的谷地线提取算法。
该算法的思想是，为每个网格赋予单元量的水，再让水按照最陡坡下流直至边界，同时将该单元量累加到流经的网格上，接着给定阈值，将累积量大于该阈值的网格视为谷地线。
洼地的处理是流水累积法的关键步骤之一，因为洼地的存在会使得某些点的径流路径发生中断，进而造成洼地点上游流水无法传递至下游。
该环节上，本文遵循以往的研究，假设洼地的出现是因为数据的测量误差或者表达精度不够造成，将洼地视为需要改正的数据，对它进行抬升，使之成为平地，以待进一步处理中赋予流向。
图2给出了在垂直地面剖面上洼地点的基本形态，从中可以看出，以某点为起始点进行邻域搜索即可检测出洼地点所在洼地的范围，其中关键在于如何控制二维平面上搜索的方向。
为了解决该问题，首先需要检测出哪些点是洼地点，洼地点是无流向点的子集，它们的特点是存在比本身高的邻居但是不存在比本身低的邻居，如图2（a）、（b）和（c）中的黑色圆点，而其中的黑色三角形点也是没有流向的点，并且也有可能需要抬高，但是不作为起始洼地点，因为可以通过设定搜索策略来计入它们；第二步就是使用给定的洼地点搜索其所在洼地范围，该过程类似于向洼地中逐渐注水，当水位到达洼地边缘时停止，并将该范围高度重新设定为水位高度。
（1）给定一个足够大的数值Ef，将下一未被标识的洼地点P放入栈，并且标识之；
（2）从栈中取栈顶点Q，其高程为EQ，若EQ大于Ef，则转至第4步；否则，考察该点每一个未被标识的邻居点M，其高程为EM，若EM小于Ef且EM不低于EQ，则点M入栈，并标识之，若EM小于Ef且EM小于EQ，则令Ef＝EQ；
（3）选择当前栈中高程最低的点，将其移动到栈顶，返回至第2步；
（4）考察所有经历过进栈和出栈的点，如果其高程低于Ef，则令其高程为Ef，点P处理完毕，其所在的洼地范围即为高程发生变化的所有点的区域。
在对所有的洼地点做上述处理之后，所有洼地都被填平。
这里没有区分以往研究中提及的独立洼地、嵌套洼地、复合洼地等情况，在洼地点搜索其洼地范围过程中，这些类型的洼地都可以被检测并处理。
对于图2（e）中的情形，如果先处理洼地点C，可以确定其范围为ABCDE，其中B、C和D被抬升至E的高度，当处理点G时，其洼地范围是ABCDEFGH，所有低于H的点都被抬升。
因此这是一个逐渐扩充的处理算法，邻接洼地或者嵌套洼地的合并蕴含在搜索的过程中，同时，搜索的扩展总是从新范围内的最低点开始的，最后的检测到的范围最多比实际的洼地范围外扩一个网格，因此其效率也是很高的。
另外，如果在处理之前对所有洼地点进行排序，然后按增序的方式检测洼地范围，可以进一步提高效率。
上述过程还可以记录洼地边缘上的溢出点（可能有多个）。
H
F
D
B
A
G
E
C
（e）
（c）
（d）
（b）
（a）
图2洼地的基本形态，黑色圆点为洼地点
下一步需要确定每个栅格的流水方向，现有的方法主要包括单流向法（O’Callanghan等，1984）、无限多流向法、有限多流向（Fairfield，1991；Torboton，1997），作为最先提出的方法，单流向法（最陡坡降方向）因其简便性和适应性而广泛应用于基于流水累积提取谷地信息的研究和应用中，同时每个栅格的唯一流向有利于汇水区、谷地线拓扑关系等其他谷地信息的确定。
本文算法使用的是该方法计算网格点的流向，给定点P的流向dir（P）＝max{（EP－Ei）/Li，i∈[1，8]}，EP和Ei分别是点P和其第i个邻居的高程，（EP－Ei）/Li是点P相对于邻居点的坡度，当邻居i是P的四方向邻居时，L＝1，i是P的八方向角相邻邻居时，L＝
，如果相对于多个邻居的坡度相同，则优先选择四方向邻居，若同为四方向或者八方向邻居，则选择对应坡度大的邻居，若仍旧相同，则按照计算顺序选择第一个。
本步骤中另外一个关键问题是平地点（即最大坡度为0的点）的流向确定问题，这里使用了Martz等（Martz等，1998）的方法，也即认为不存在绝对水平的地表，并且地形表面上的流水总是沿最短路径流出平地区域，为平地点逐渐累加微量即可在平地区域形成由高到底的微小起伏。
根据流向信息和流量分配策略，从邻居入流数为0的格网点出发，将当前点上的水量添加到流向点上，将流向点作为当前点迭代该过程直至达到邻居入流数大于1或者格网边界上的点。
最后各点上的累加信息构成了流水累积量矩阵，每个累积量值表示对应点的汇水面积。
图3（b）是流水累积矩阵示例，可以看出，位于谷地的网格点所得到的累积量比局部坡面上的累积量要大，而位于高程较低区域的累积量相对的要大于较高区域。
2.2.谷地线的提取以及相关信息的构建
用给定阈值对流水累积量矩阵进行“二值化”，累积量大于阈值的点即视为谷地点，图3（b）中的数值表示对应点的累积量，当阈值设定为50时，便得到着色网格点表示的谷地点，从图中累积量数值大小的分布可以看出，设定的阈值越小，得到的谷地点越多，也就能标识出越多的特征信息。
利用流向信息可以进一步将过滤出的谷地点连接成为谷地线，以此为基础还可以得到更多有意义的信息，第一个就是谷地线的拓扑关系以及谷地线段的属性，第二个是区域内汇水区的分割信息，也即结构化信息，以往的研究都是将这些工作分步计算得到的，本文设计实现了新的算法（如图4），可以同时获取并构建谷地线和汇水区的信息。
谷地线拓扑关系的表达可以模拟河流的分级编码方法（毋河海等，1997），其中Horton编码、Strahler编码、Shreve编码等是最为常用的几种形式，本节采用了Strahler编码。
（a）高程格网及各点流向
（b）流水累积量和阈值过滤
（c）谷地线的Strahler编码
（d）汇水区划分
图3谷地线及相关信息的提取过程
算法的执行需要各点的流向信息和谷地点的邻居流入数信息，最后输出记录了谷地点Strahler编码矩阵、谷地点所在段标识矩阵，假设谷地点P的Strahler编码为SP、段标识码为OP、当前邻居入流数为IP，具体的算法流程如图4。
该算法的结果除记录了每条谷地段上游和下游顶点外，主要内容仍然是栅格矩阵，这些信息已经足够说明谷地段之间的上下游关系，为了得到谷地段的几何属性信息（例如长度、坡度等）以及汇水区划分信息，还需要根据流向信息和谷地段标识码信息进一步处理，这里只需从标识码矩阵中任意没有赋值的点开始，沿流向前进直至有标识码的点，则两点之间所有点都被赋值为该标识码，最后得到的矩阵就是汇水区矩阵，其中每个点都被赋予其直接流入的谷地段标识码，该值即为对应谷地段的汇水区标识，同时可以统计各个汇水区的面积，应当注意的是，图4算法中谷地段的标识码是从1开始的，因此汇水区矩阵中0值对应的就是这些流向区域外的部分。
对于谷地段的长度，可以使用Strahler编码矩阵和汇水区矩阵来统计得到，而坡度可以从记录的上下游结点和长度计算出来。
2.3.局部形态法扩展基本谷地线
流水累积法的全局性克服了固定窗口只能检测有限范围的缺点，但是阈值分割使得该方法会丢失较高区域的谷地线，其原因阈值无法在局部自适应的调整分割，这样该方法又显得过于“全局性”，在图像处理、人工智能等研究中的可变阈值方法在这里并不适用，因为，一则各网格点上流水累积量的数值范围与图像中的灰度范围有较大的差异，二则是谷地系统有着自身独特的自然特性，这些可变阈值或者自适应阈值的确定并没有顾及到这一点。
局部形态法由于只是考虑较小范围内的相对程度，因此可以弥补流水累计法的不足。
局部形态法扩展过程分为不完整谷地线上游顶点延伸和形态谷地点下行两个步骤。
首先，流水累积法提取的某些谷地段相对于实际情况而言是不完整的，在其上游顶点之上还存在着形态十分明显的谷地。
以往基于形态法提取谷地线的研究中，大多是采用从上游点向下游连接的方式（Seemuller，1989；Takahashi等，1995），本文则是以已有的部分谷地线为启发性知识，由其上游顶点向上连接合适网格点以补全丢失的部分，具体的实现过程中，首先取一条一级谷地段的上游顶点，然后搜索满足流向当前点、是局部形态的谷地点、当前不是谷地点等三个条件的邻居点，将其标识为谷地点，接着再以这些新增的谷地点为新的出发点继续做同样检测，直至找不到满足条件的邻居点。
N
N
Y
N
N
若OP>0，则OID＝OP
OP＝OID
N
Y
Y
SP＝SID
由流向信息更新P为下一点
SID＝SP，记录P
SID＝SID+1，记录P;StrhSta[SID]++
OID＝OID＋1；OP＝OID
初始化整型变量SID=1;StrhSta[SID]++
Y
Y
N
OID＋＋；OP＝OID，记录P
SP＝SID
SP＝SID+1
N
退出
初始化整型变量OID＝1
Y
Y
图4构建谷地线的Strahler编码和标识码:
给定点若有N个谷地点邻居流入，则其邻居入流数记为N＋1
其次，上步处理仅补齐了不完整谷地线，而且需要已有的谷地段作为引导，然而，还有大量在流水累积法的结果没有任何体现的谷地。
不过我们注意到当前的结果已经基本上包含了谷地线主干，未检测出的谷地段恰是形态上比较典型而因为处于较高区域，因此完全可以采用由上而下的连接策略，本步处理实现中的基本思想是，以当前不属于谷地而局部形态法检测出是谷地的点为起点，沿最陡坡向下行，直至已有的谷地点。
为了得到正确合理的连接结果，需要考虑多个因素。
流水累积法的结果已经包括了宽阔、低缓区域的谷地线，因此这些区域的谷地基本上无需再次补充，同时形态法也不适应于起伏缓和区域，因此对下行连接起点的选择施加一定的坡度限制，也即坡度小于给定值的点不能作为起点；由于下行连接过程中所经各点地形特征在整体表现上会存在很多情况，这里采用三种手段降低可能出现的不确定性。
首先，网格点的邻居入流数能说明其作为某种地形特征的程度，入流数较大的点多分布在谷地形态显著的位置，所以一个谷地点的邻居入流数越大，则连接成为谷地线的可能性越大。
其次，每点下行中的方向整体上是沿着最陡坡向，然而如果该方向指的是山脊点，则重新在邻居依次搜索比当前点低的局部形态谷地点或者坡面点，若存在，该修改流向指向新点；如果最陡坡向指向坡面点，做类似处理。
第三条策略是，在下行过程中分别统计其中谷地点和山脊点的数目，下行结束后，将它们分别作为有效长度和无效长度，只有当有效长度大于无效长度、有效长度与全长比值大于给定值、无效长度与全长比值小于给定值时，才真正将下行路径上各点标记为谷地点。
3.谷地线的形态整理与噪声消除
如前文所述，流水累积法提取的谷地线在形态和关系上存在需要改进的地方，实际上该方法的应用过程中，洼地和平地的存在、结果中谷地线的粘连以及平直谷地线都是得到合理结果的障碍。
而使用形态法扩展该结果会引入更多问题，主要体现为将不符合谷地点特征的点归为谷地点、将符合其特征的点漏选以及谷地线之间的粘连。
本文将这些问题通称为噪声，归纳而言，洼地和平地的无流向、符合特征点的漏选等问题已经基本解决，下面将首先分析平直谷地线和粘连谷地线的特点，对它们进行分类，进而设计相应的规则和对策修正。
这里首先定义粘连，在谷地段矩阵中，若有两个谷地点流向互不指向，当它们是四方向相邻时，则称它们为粘连，当它们是八方向对角相邻时，则称它们为轻度粘连。
从谷地段级别的角度看，它们之间存在着同级之间、不同级之间的粘连，而以谷地段粘连部分长度为标准，分为全段粘连、下游部分粘连、上游部分粘连、中间部分粘连、环状粘连型等。
对这些粘连进行清理校正过程中，主要假设是认为在Strahler编码的谷地段中，级别低的重要性低于级别高的，其次参考流水累积量、长度等谷地段属性。
3.1.平直谷地线
平直谷地线的出现使得提取结果显得非常不自然，同时它们改变了区域谷地密度的对比，分析平直谷地线的特点，可以看出它们出现的基本特点：
（1）主要出现在地形变化很平缓、坡度很小的区域，比如宽阔的谷地；
（2）平直谷地线本身的坡度很小；（3）谷地线上各个点的流向基本上是相同的；（4）谷地段上形态谷地点的数目比例很小；（5）谷地线从上游到下游的各点流量累积变化值比较固定而且很小，也即相邻的两个谷地线点上的累积流量相差很小；（6）谷地线的Strahler级别大多为1，当流水累积分割值适当的时候，它们的长度都是比较短的。
在分析了平直谷地线的特点之后，可以以此为标准对应的删除这些点，不过我们应该看到，对于其中一些特点的判定是难以给出适应所有情况的参数，比如坡度。
3.2.长度为1的一级谷地段
经过局部形态法扩展后，结果中存在大量的长度为1的一级谷地段，而其中又有相当一部分是与其他谷地段粘连在一起。
该处理将满足下列条件之一的视为噪声并删除：
（1）四方向上存在更高级别的邻居，且该邻居不是当前点的下游点；
（2）四方向上存在两个或者以上的邻居，且都不是当前点的下游点。
图5是该处理前后的对比效果。
3.3.粘连于高级别谷地段的一级谷地段
这种情况是一级谷地段上若干点粘连与其他级别谷地段上，从一级谷地段的角度出发，可以分为下游粘连和上游粘连。
（a）原始图
（b）删除后
图5长度为1的噪声谷地段（红色为1级谷地）
图6一级谷地段粘连与高级谷地，图中数字为对应点流水累积量
下游粘连是从一级谷地段下游顶点开始就粘连（包括轻度粘连）与高级别谷地点上，对于这种情况则从下游点开始上溯逐点删除，直至当前点不再与高级谷地段粘连，这里会有两种结果，一是当前段上所有点都发生粘连，那么都被删除，如图6左图下部矩形框内的情形；第二是当前段上游某点不粘连，图6左图上部矩形框内即是此情形，这时只有下游部分谷地点被删除（即图中标记11的网格点），为了使得剩余的部分不至于丢失拓扑关系，对最后一个粘连点做特殊处理，也即为它寻找一个合适的下游邻居点，这个下游点应当满足的条件：
低于当前点、是一个谷地点、是当前点邻居中汇流量最大的，图6左图上部矩形内粘连的最后点是标记6的点，在点11删除后该点将指向标记35的点，同时在6原来指向的路径上删除本点的累积量（6），而在新指向的路径上增加该值。
上游粘连是除了下游顶点外的上游某点与其他级别谷地段发生粘连（非轻度粘连），对于这种情况，本文将该段在粘连处分解成两段，图6右图矩形内即为这种情况，该点将根据上一段中修改流向的方法重新寻找一个下游点，并且调整新旧路径上的流水累积量。
3.4.一级谷地线段下游相互粘连
这种情形还需要细分为三种情况，第一种是两根一级谷地线段相互粘连，其具体特征是两条一级谷地段的下游顶点是相互粘连的，而且粘连向上游方向持续的长度大于1，如果只有下游顶点粘连，则不视为谷地段粘连。
为了消除粘连，需要选择其中之一删除，标准是优先保留流水累积量大、长度大、有效长度大、坡度大、下游点流向与下游段方向的相近程度大的谷地段，对五个指标依次判断，如果都相同则任意删除其中之一。
这里同样会出现一条谷地段部分被删除的情形，对剩余部分的处理与3.3节中的处理相同。
第二种情况是三根一级谷地线段在下游相互粘连，本文在处理时没有深入考虑汇水量等特征，而是删除两侧、保留中间，若被删除的谷地存在剩余段，则修正其最后一个粘连点。
第三种是多于三根的一级谷地线段在下游相互粘连，类似于第二种情况的处理，这里仍旧没有考虑相关的物理特征，仅仅是为了得到单象素宽的结果，采用了尽可能简单的方法，具体处理过程中，锁定这些相互粘连区域两侧的谷地段，这两条谷地段都只有一个邻居谷地段，然后删除紧邻的两条谷地段的粘连。
对剩余部分做同样处理。
图7两根一级谷地线段粘连
图8两根以上一级谷地段下游粘连
3.5.两根一级谷地段上游粘连
图9高级别谷地段之间的粘连
升级会员 