CN107368839B - 一种基于dem的断层面的自动提取方法 - Google Patents

一种基于dem的断层面的自动提取方法 Download PDF

Info

Publication number
CN107368839B
CN107368839B CN201710479465.3A CN201710479465A CN107368839B CN 107368839 B CN107368839 B CN 107368839B CN 201710479465 A CN201710479465 A CN 201710479465A CN 107368839 B CN107368839 B CN 107368839B
Authority
CN
China
Prior art keywords
sequence
polygon
variance
dem
file
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
CN201710479465.3A
Other languages
English (en)
Other versions
CN107368839A (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.)
Nanjing Normal University
Original Assignee
Nanjing Normal 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 Nanjing Normal University filed Critical Nanjing Normal University
Priority to CN201710479465.3A priority Critical patent/CN107368839B/zh
Publication of CN107368839A publication Critical patent/CN107368839A/zh
Application granted granted Critical
Publication of CN107368839B publication Critical patent/CN107368839B/zh
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
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/40—Extraction of image or video features
    • G06V10/46—Descriptors for shape, contour or point-related descriptors, e.g. scale invariant feature transform [SIFT] or bags of words [BoW]; Salient regional features
    • G06V10/462—Salient features, e.g. scale invariant feature transforms [SIFT]

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Multimedia (AREA)
  • Theoretical Computer Science (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于DEM的断层面的自动提取方法。包括以下几个步骤：S1、根据输入的DEM数据，进行坡度矩阵生成；S2、根据生成的坡度矩阵，进行初步区域筛选；S3、根据初步筛选出的区域，计算其区域方差，进行二次区域筛选。本发明利用了DEM数据，直接根据研究所得的断层面的一些性质进行自动化提取，因此相比较目前基于野外调查或遥感影像，人工化进行识别和提取的方法，自动化程度和效率都比较高。

Description

一种基于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′；
Figure GDA0001412365240000021
(其中i，j为元素的行号和列号，a1,a2,a3,a4,a5,a6,a7,a8分别为该元素的左上方、上方，右上方、左方、右方、左下方、下方、右下方对应元素的值，sij为该点的坡度值)；
Q2、根据生成的坡度矩阵A′，进行初步区域筛选；
Q21：坡度矩阵降维排序：将坡度矩阵A′以行序为主序,变换为一维序列B；再将一维序列B降序排序得到一维序列B′；
Q22：令
Figure GDA0001412365240000022
表示序列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)，
Figure GDA0001412365240000041
(其中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数据和断层面本身存在的数理性质，因此本方法具有经济性、高效性和自动化程度高等特点。

Claims (1)

1.一种基于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′；
Figure FDA0002406098880000011
其中i，j为元素的行号和列号，a1,a2,a3,a4,a5,a6,a7,a8分别为该元素的左上方、上方，右上方、左方、右方、左下方、下方、右下方对应元素的值，sij为该元素中心点的坡度值，b表示坡度值在二维矩阵列方向的分量，c表示坡度值在二维矩阵行方向的分量；cellsize为像元大小；
Q2、根据生成的坡度矩阵A′，进行初步区域筛选；
Q21：坡度矩阵降维排序：将坡度矩阵A′以行序为主序,变换为一维序列B；再将一维序列B降序排序得到一维序列B′；
Q22：令
Figure FDA0002406098880000012
表示序列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所包含的多边形，即为提取出的断层面。
CN201710479465.3A 2017-06-22 2017-06-22 一种基于dem的断层面的自动提取方法 Active CN107368839B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710479465.3A CN107368839B (zh) 2017-06-22 2017-06-22 一种基于dem的断层面的自动提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710479465.3A CN107368839B (zh) 2017-06-22 2017-06-22 一种基于dem的断层面的自动提取方法

Publications (2)

Publication Number Publication Date
CN107368839A CN107368839A (zh) 2017-11-21
CN107368839B true CN107368839B (zh) 2020-05-12

Family

ID=60304963

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710479465.3A Active CN107368839B (zh) 2017-06-22 2017-06-22 一种基于dem的断层面的自动提取方法

Country Status (1)

Country Link
CN (1) CN107368839B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110413712B (zh) * 2019-06-11 2021-09-28 南京泛在地理信息产业研究院有限公司 一种基于dem的断层面自动识别方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101936008A (zh) * 2010-09-30 2011-01-05 东北大学 岩体边坡三维模型及块体滑落分析方法
CN102819023A (zh) * 2012-07-27 2012-12-12 中国地质大学（武汉） 基于LiDAR的复杂地质背景区滑坡识别的方法及系统
CN103412331A (zh) * 2013-08-30 2013-11-27 电子科技大学 一种三维地震断层自动提取方法
CN104166163A (zh) * 2014-08-27 2014-11-26 电子科技大学 基于三维大数据量地震数据体的断层曲面自动提取方法
CN106649466A (zh) * 2016-09-27 2017-05-10 西安电子科技大学 数字地图中典型地形几何参数获取方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9618639B2 (en) * 2012-03-01 2017-04-11 Drilling Info, Inc. Method and system for image-guided fault extraction from a fault-enhanced seismic image

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101936008A (zh) * 2010-09-30 2011-01-05 东北大学 岩体边坡三维模型及块体滑落分析方法
CN102819023A (zh) * 2012-07-27 2012-12-12 中国地质大学（武汉） 基于LiDAR的复杂地质背景区滑坡识别的方法及系统
CN103412331A (zh) * 2013-08-30 2013-11-27 电子科技大学 一种三维地震断层自动提取方法
CN104166163A (zh) * 2014-08-27 2014-11-26 电子科技大学 基于三维大数据量地震数据体的断层曲面自动提取方法
CN106649466A (zh) * 2016-09-27 2017-05-10 西安电子科技大学 数字地图中典型地形几何参数获取方法

Also Published As

Publication number Publication date
CN107368839A (zh) 2017-11-21

Similar Documents

Publication Publication Date Title
CN106815847B (zh) 基于激光雷达点云的树木分割方法及单棵树提取方法
Nie et al. A revised progressive TIN densification for filtering airborne LiDAR data
Ke et al. A review of methods for automatic individual tree-crown detection and delineation from passive remote sensing
Guan et al. Integration of orthoimagery and lidar data for object-based urban thematic mapping using random forests
CN104134234B (zh) 一种全自动的基于单幅图像的三维场景构建方法
CN112396128B (zh) 一种铁路外部环境风险源样本自动标注方法
Li et al. A spatial–temporal Hopfield neural network approach for super-resolution land cover mapping with multi-temporal different resolution remotely sensed images
CN102930518B (zh) 基于改进的稀疏表示的图像超分辨率方法
CN111462838B (zh) 一种直接将图像像素转换成有限元单元的方法
Iosifescu et al. Towards a comprehensive methodology for automatic vectorization of raster historical maps
CN108986024A (zh) 一种基于网格的激光点云规则排列处理方法
Rosim et al. TerraHidro: a distributed hydrology modelling system with high quality drainage extraction
Rashidi et al. Ground filtering LiDAR data based on multi-scale analysis of height difference threshold
CN107368839B (zh) 一种基于dem的断层面的自动提取方法
CN104992407A (zh) 一种图像超分辨方法
CN114842356B (zh) 一种高分辨率地表类型样本自动生成方法、系统及设备
WO2024020744A1 (zh) 基于高分辨率影像的生态图斑提取方法
Erfanifard et al. A robust approach to generate canopy cover maps using UltraCam-D derived orthoimagery classified by support vector machines in Zagros woodlands, West Iran
Ai et al. Improved sub-pixel mapping method coupling spatial dependence with directivity and connectivity
CN115410036A (zh) 一种高压架空输电线路关键要素激光点云自动分类方法
CN107220651A (zh) 一种提取图像特征的方法及装置
CN110288645B (zh) 一种基于坡向约束的地形表面面积计算方法
CN107958485B (zh) 一种伪山顶点剔除方法
CN113421194A (zh) 一种根据布格重力异常图像提取隐伏断层的方法
Hassan Dynamic expansion and urbanization of greater Cairo metropolis, Egypt

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