一种基于DEM的断层面的自动提取方法与流程

本发明属于地理信息自动解析领域，具体涉及一种利用dem自动提取断层面的方法。

背景技术：

断层面是指断层崖经河流或冲沟切割侵蚀后形成的三角形陡崖，是现代活动断层的标志。断层面自动识别和提取，对于地质资源勘探、工程建设以及地质环境评价等工作具有重要意义。

dem作为一种基础地理信息资源，能够较好表现原生地貌特征。长期以来，基于dem已能够有效进行坡度、坡向、山体阴影等信息的自动化提取。断层面作为一种常见的地貌类型，目前，还主要是基于野外调查或遥感影像，人工化进行识别和提取，相关工作投入大、效率较低，且解析质量难以保证。为此，迫切需要研发一种经济、高效的断层面自动化提取方法。

技术实现要素：

本发明的目的在于，利用dem提供的地形特征和断层面的几何特征，通过坡度矩阵生成、基于坡度的区域筛选和基于方差的断层面提取三个步骤，提供一种自动化提取断层面的方法。

为了实现上述目的，本发明采取的技术方案如下：

一种基于dem的断层面的自动提取方法，包括以下几个步骤：q1、根据输入的dem数据，生成坡度矩阵a′；

q11：dem数据预处理：读取出dem的行数m，列数n，像元大小与空值默认值nodata，进一步将dem数据按照行序为主序，转变为二维矩阵a；

q12：提取参数设置：根据dem的数据信息设置相应的提取参数，提取参数包括输入排序百分比参数σ％(σ∈(1,100))、面聚合距离distance、面积选择百分比参数ρ％(ρ∈(1,100))和方差因子t；

q13：坡度计算：判断二维矩阵a中的每一元素的左上方、上方，右上方、左方、右方、左下方、下方、右下方是否存在元素，且都不为空值默认值nodata；若成立，则根据公式(1)计算出的坡度矩阵a′；

(其中i，j为元素的行号和列号，a1,a2,a3,a4,a5,a6,a7,a8分别为该元素的左上方、上方，右上方、左方、右方、左下方、下方、右下方对应元素的值，sij为该点的坡度值)；

q2、根据生成的坡度矩阵a′，进行初步区域筛选；

q21：坡度矩阵降维排序：将坡度矩阵a′以行序为主序,变换为一维序列b；再将一维序列b降序排序得到一维序列b′；

q22：令表示序列b′中的下标；遍历a′，若a′ij≥b′v，则将其值设置为1，否则设置为0，形成新的坡度矩阵a″；

q23：将新的坡度矩阵a″变换为栅格数据后转换为矢量文件polygon；

q24：面聚合与填岛:根据多边形面聚合算法，对矢量文件polygon中的距离小于面聚合距离distance的面进行聚合操作，得到矢量文件polygon′；利用多边形左转算法，对矢量文件polygon′进行填岛操作，得到新的矢量文件polygon″；

q25：按面积筛选多边形：将新的矢量文件polygon″中包含的多边形按照面积降序排列，并抽取前ρ％个多边形，组成新的多边形序列t，t＝{p1,p2,…，pk}。

q3、根据初步筛选出的区域，计算其区域方差，进行二次区域筛选；

q31：将坡度矩阵a′转换为栅格文件slope；

q32：区域方差计算：对多边形序列t中的每一个多边形pi，循环执行1)-3)的步骤，得到区域方差序列v；

1)利用pi对slope进行掩膜运算，得到栅格文件mask；

2)以行序为主序读取栅格文件mask，提取出非空的值，形成栅格文件mask的栅格值序列z＝{z1,z2,…,zn}(其中n为栅格文件mask包含的栅格单元的数量)；

3)计算序列z的方差vi，得到元组(pi,vi),并插入区域方差序列v中；

q33:对区域方差序列v中的每一元组，若vi＜t(t为方差因子),则将元组(pi,vi)插入到新的方差序列v′中；

q34：遍历方差序列v′，读出其中的所有多边形集合，并写入一个新的矢量文件shp，该新的矢量文件shp所包含的多边形，即为提取出的断层面。

有益效果：本发明利用了dem数据，直接根据研究所得的断层面的一些性质进行自动化提取，因此相比较目前基于野外调查或遥感影像，人工化进行识别和提取的方法，自动化程度和效率都比较高。

附图说明

图1为本发明的流程图。

图2为实验数据。

图3坡度筛选后的栅格数据。

图4栅格转矢量后的文件。

图5面聚合后的文件。

图6填岛后的文件。

图7多边形序列t。

图8mask文件。

图9本方法提取结果与实测断层线的对比图。

具体实施方式

下面结合本发明的附图，对本发明的技术方案进行清楚、完整地描述。显然所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1-9所示，本发明实施例提供一种基于dem的断层面提取方法，如下：

s1、根据输入的dem数据，进行坡度矩阵生成，具体包括如下步骤：

(一)dem数据预处理

将本实施例选取的dem数据在arcscene中显示如图2所示。

基于arcengine软件，读取出dem的行数m＝935，列数n＝1220，像元大小cellsize＝5与空值默认值nodata＝-9999，进一步将dem数据按照行序为主序写入本实例开辟的935行，1220列的二维矩阵a中；

(二)提取参数设置

输入排序百分比参数15％。输入面聚合距离50m。输入面积选择百分比参数为15％。输入方差因子为25；

(三)坡度计算

以行数m＝50,列数n＝150，高程值为831的点为例。其左上方、上方，右上方、左方、右方、左下方、下方、右下方的元素都存在，且值分别为835,838,834,837,827,829,818,815；根据公式(1)，

(其中i，j为元素的行号和列号，a1,a2,a3,a4,a5,a6,a7,a8分别为该元素的左上方、上方，右上方、左方、右方、左下方、下方、右下方对应元素的值，sij为该点的坡度值)。

得到b＝-0.375，c＝1.625，s＝1.0307；对矩阵a中的每一元素都进行如上操作，得到新矩阵a′；

s2、根据生成的坡度矩阵，进行初步区域筛选

(一)坡度矩阵降维排序

将二维矩阵a′以行序为主序,变换为一维序列b＝{-9999,-9999,…,-9999}；再将b降序排序得到一维序列b′＝{86.5,86.4…,-9999}；

(二)计算得到v＝171105，a″ij＝{0,0,0,1,…，0}；

(三)基于arcengine软件，将a″变换为栅格数据，在arcmap软件中显示如图3，将此栅格数据转换为矢量文件polygon，在arcmap软件中显示如图4。

(四)面聚合与填岛

根据多边形面聚合算法，对矢量文件polygon中的距离小于50m的面进行聚合操作，得到矢量文件polygon′，在arcmap中显示如图5；

利用多边形左转算法，对矢量文件polygon′进行填岛操作，得到polygon″，在arcmap中显示如图6；

(五)按面积筛选多边形

将矢量文件polygon″中包含的多边形按照面积降序排列，并抽取前15％个多边形，组成新的多边形序列为t，在arcmap中显示如图7；

s3、根据筛选出的区域，计算其区域方差，进行二次区域筛选；

(一)基于arcengine软件将坡度矩阵a′转换为栅格文件slope；

(二)区域方差计算：

1)以t中第19个多边形为例，基于arcengine，对slope进行掩膜运算，得到栅格文件mask，栅格文件mask文件在arcmap中显示如图8；

2)以行序为主序读取栅格文件mask文件，提取出非空的值，形成栅格文件mask的栅格值序列z＝{71.6,71.6,…,70.5}(其中n为mask包含的栅格单元的数量)；

3)计算序列z的方差，得到元组(p19,21.921),并插入区域方差序列v中。

4)对t中的每一元素执行步骤3)操作，得到

v＝{(p1,22.225),(p2,22.225),…,(p22,21.924)}；

(三)对v中的每一元组计算，若vi＜t,则将(pi,vi)元组插入到新的方差序列v′中，得到v′＝{(p1,22.225),(p2,22.225),…，(p19,22.001)；

(四)遍历v′，读出其中的所有多边形集合，并写入一个新的矢量文件shp。该新的矢量文件shp所包含的多边形，即为提取出的断层面，断层面在arcscene中显示如图9。

测试分析：将提取出的断层面和实测断层线进行叠置，如图9，可以看出两者具有较高的一致性，说明本方法的有效性。

此外，与传统基于遥感影像的人工解译方法相比，由于使用了dem数据和断层面本身存在的数理性质，因此本方法具有经济性、高效性和自动化程度高等特点。