CN107255516B - 一种遥感影像滑坡单体划分方法 - Google Patents
一种遥感影像滑坡单体划分方法 Download PDFInfo
- Publication number
- CN107255516B CN107255516B CN201710393741.4A CN201710393741A CN107255516B CN 107255516 B CN107255516 B CN 107255516B CN 201710393741 A CN201710393741 A CN 201710393741A CN 107255516 B CN107255516 B CN 107255516B
- Authority
- CN
- China
- Prior art keywords
- landslide
- monomer
- common boundary
- remote sensing
- threshold
- 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
Links
- 239000000178 monomer Substances 0.000 title claims abstract description 87
- 238000000034 method Methods 0.000 title claims abstract description 84
- 238000000605 extraction Methods 0.000 claims abstract description 47
- 230000004927 fusion Effects 0.000 claims abstract description 30
- 238000002156 mixing Methods 0.000 claims abstract description 5
- 239000000284 extract Substances 0.000 claims description 10
- 230000011218 segmentation Effects 0.000 claims description 9
- 238000001228 spectrum Methods 0.000 claims description 7
- 230000004069 differentiation Effects 0.000 claims 1
- 238000011156 evaluation Methods 0.000 abstract description 10
- 230000003595 spectral effect Effects 0.000 description 11
- 230000000007 visual effect Effects 0.000 description 11
- BVPWJMCABCPUQY-UHFFFAOYSA-N 4-amino-5-chloro-2-methoxy-N-[1-(phenylmethyl)-4-piperidinyl]benzamide Chemical compound COC1=CC(N)=C(Cl)C=C1C(=O)NC1CCN(CC=2C=CC=CC=2)CC1 BVPWJMCABCPUQY-UHFFFAOYSA-N 0.000 description 6
- 238000012545 processing Methods 0.000 description 6
- 238000012937 correction Methods 0.000 description 4
- 238000012795 verification Methods 0.000 description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 4
- 230000006399 behavior Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000003709 image segmentation Methods 0.000 description 3
- 239000000126 substance Substances 0.000 description 3
- 238000013461 design Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 239000011435 rock Substances 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 239000004576 sand Substances 0.000 description 2
- 238000007619 statistical method Methods 0.000 description 2
- 241000208340 Araliaceae Species 0.000 description 1
- 241000196324 Embryophyta Species 0.000 description 1
- 240000007594 Oryza sativa Species 0.000 description 1
- 235000007164 Oryza sativa Nutrition 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000009313 farming Methods 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- VMXUWOKSQNHOCA-UKTHLTGXSA-N ranitidine Chemical compound [O-][N+](=O)\C=C(/NC)NCCSCC1=CC=C(CN(C)C)O1 VMXUWOKSQNHOCA-UKTHLTGXSA-N 0.000 description 1
- 238000002310 reflectometry Methods 0.000 description 1
- 235000009566 rice Nutrition 0.000 description 1
- 230000001932 seasonal effect Effects 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J3/00—Spectrometry; Spectrophotometry; Monochromators; Measuring colours
- G01J3/28—Investigating the spectrum
- G01J3/2823—Imaging spectrometer
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J3/00—Spectrometry; Spectrophotometry; Monochromators; Measuring colours
- G01J3/28—Investigating the spectrum
- G01J3/2823—Imaging spectrometer
- G01J2003/2826—Multispectral imaging, e.g. filter imaging
Landscapes
- Physics & Mathematics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Theoretical Computer Science (AREA)
- Image Analysis (AREA)
Abstract
本发明提供一种遥感影像滑坡单体划分方法。该方法包括：以研究区域的多光谱影像数据和DEM数据为基础，进行滑坡提取得到滑坡对象；设置NDVI标准差阈值，利用融合算法将所得到的滑坡对象融合；利用滑坡间的共同边界长度阈值使经融合的滑坡对象进一步融合得到滑坡主干；以及利用共同边界占滑坡周长百分比阈值使得滑坡主干与其周围的滑坡对象融合，得到滑坡单体。根据本发明的滑坡单体划分方法则可以将提取得到的滑坡对象划分成一个个的滑坡单体，避免出现不合理的巨大滑坡区域。特别地，可对地震或强降雨引起的大范围滑坡灾害进行快速评估，得到滑坡灾害的面积与数量，为灾后应急救援和恢复重建工作提供参考。
Description
技术领域
本发明涉及遥感影像数据处理。更具体地，本发明涉及一种基于遥感影像的滑坡单体划分方法。
背景技术
滑坡是组成斜坡的岩石、土壤、人工填埋物或多种物质的向下和向外移动，包括典型的滑动、落石和流动。地震是诱发滑坡的主要因素，并常常造成大量的人员伤亡和经济损失。2008年的汶川地震诱发了大量滑坡等地质灾害，不仅破坏了道路、输电线、管道等生命线工程，掩埋了大量农田与森林，还造成了约20000名人员伤亡。因此，无论是从震后灾情的迅速评估、应急救援力量的合理配置，还是灾后恢复重建工程的科学选址方面考虑，都需要在灾后迅速掌握滑坡灾害的情况。与传统的人工排查大范围滑坡灾害相比，通过遥感影像提取大范围滑坡具有危险性小、速度快、少遗漏、覆盖范围广、可实现半自动或全自动计算机处理等特点，因此得到广泛应用。
地震以及强降雨引发的滑坡往往成片分布，众多的滑坡体之间会通过小的连通区域形成结构复杂，面积巨大的滑坡区域。例如汶川地震形成了大量彼此相连，结构复杂，面积巨大的滑坡区域，给滑坡单体的识别增大了困难。计算机自动提取滑坡具有速度快，效率高，不受主观因素影响的特点，被越来越多的学者所青睐。滑坡的面积与数量是反映滑坡灾害严重程度的重要指标，面积反映滑坡的范围大小，数量则可以表现滑坡的密集程度。要准确统计滑坡个数就必须对这些面积巨大的滑坡区域进行单体划分。
随着遥感影像提取滑坡的广泛应用，滑坡单体划分将成为这一领域的一个难题。因此，需要提供一种可应用于较大范围地区遥感影像滑坡提取方法与滑坡对象的单体划分方法，为灾后快速而准确地获得大范围滑坡灾害的面积与数量提供技术支持。
发明内容
为解决上述问题，本发明提出一种基于遥感影像的滑坡单体划分方法。彼此相连的不同滑坡，由于空间上紧挨在一起，其滑坡碎屑物相似性较大。研究区域的滑坡绝大多数是震后几个月内发生的新滑坡，由于发生的时间相近，滑坡面有很多相似的地方，这为滑坡的单体划分增加了难度。本发明认为滑坡体内部物质间的光谱差异较小，而滑坡体之间的连通区域有较大的光谱差异，突出表现在NDVI上。滑坡的连通区域位于滑坡体的边沿，有限的碎屑物不能把植被完全掩埋，因此，与滑坡主干内部相比，连通区域存在较多NDVI较高的象元。对此，本发明在融合滑坡对象时，通过设置NDVI的标准差阈值，使得滑坡主干与连通区域各自融合成较大的对象。在几何特征上，滑坡连通区域往往较窄，滑坡间的共同边界长度占滑坡体周长的比例较低。利用这一特征，在滑坡单体划分时，通过设置滑坡对象间共同边界的长度阈值、共同边界长度占滑坡周长百分比阈值将不同滑坡体分开。
根据本发明的遥感影像滑坡单体划分方法包括：
以研究区域的多光谱影像数据和数字高度模型(DEM，Digital Elevation Model)数据为基础，进行滑坡提取，得到滑坡对象；
设置归一化植被指数(Normalized Difference Vegetation Index，NDVI)标准差阈值，利用融合算法将所得到的滑坡对象融合；
利用滑坡间的共同边界长度阈值使经融合的滑坡体内部滑坡对象进一步融合得到滑坡主干；以及
利用滑坡间的共同边界占滑坡周长百分比阈值使得滑坡主干与其周围的滑坡对象融合，由此得到滑坡单体。
优选地，该方法进一步包括根据滑坡面积大小将研究区域的滑坡对象分为多个等级，分别对各等级的滑坡对象设置NDVI标准差阈值、共同边界长度阈值和共同边界占滑坡周长百分比阈值，进行融合。
优选地，所述设置NDVI标准差阈值包括通过试错法设置NDVI标准差阈值，设置要求是使滑坡体滑坡对象尽可能融合，而连通区内滑坡对象不与滑坡体合并。
优选地，按滑坡对象共同边界占周长的平均百分比确定其共同边界长度阈值。
优选地，滑坡对象共同边界占周长的百分比阈值大于该滑坡对象共同边界占周长的平均百分比。
优选地，所述以研究区域的多光谱影像数据和DEM数据为基础，进行基于滑坡提取得到滑坡对象的步骤进一步包括：
对多光谱影像数据进行多尺度分割；
对经分割的影像进行滑坡候选对象的提取，
构造基于光谱、纹理、地形和上下文关系特征的规则，逐步剔除非滑坡对象，获得滑坡对象。
优选地，所述对经分割的影像进行滑坡候选对象的提取的步骤包括利用NDVI滑坡提取阈值区分植被对象和滑坡对象，提取滑坡候选对象。
优选地，所述逐步剔除非滑坡对象的步骤进一步包括，按照积雪、云、阴影、河、湖泊、人造物的顺序剔除非滑坡对象。
优选地，该方法进一步包括对得到的滑坡对象进行棋盘分割，以剔除小的植被对象。
优选地，该方法进一步包括，利用Find Enclosed by Class算法，将被滑坡体包围的滑坡对象融合在一起。
本发明在较大范围研究区域内构建了基于遥感影像的滑坡提取方法。首先设定NDVI阈值将植被对象剔除得到滑坡候选对象。随后运用各种方法，根据影像的光谱、纹理、地形等特征构建土壤-水体指数(LWM)、灰度共生矩阵(GLCM)、横向系数、长宽比等指标，将常年积雪区、季节性积雪区、云、阴影、河流、河沙、湖泊、道路、建成区、农用地这些非滑坡对象逐类剔除。之后对提取得到的滑坡对象进行小尺度的棋盘分割，进一步剔除植被对象，得到最终的滑坡对象并融合。
本发明进一步提出了滑坡单体划分方法。通过对剔除植被象元后的棋盘分割对象进行多尺度分割融合，得到与棋盘分割前形状类似的滑坡对象。根据滑坡内部NDVI差异小，而滑坡之间连接区NDVI差异大的特点，通过构建NDVI标准差阈值使滑坡内部融合成较大的对象。通过将研究区域的滑坡体按面积从大到小分成多个等级，设定各滑坡对象共同边界长度阈值使滑坡内部较大的对象合并成滑坡主干，再设定共同边界占滑坡周长百分比阈值，使滑坡主干与其周边的对象合并形成完整的滑坡单体。
根据三个验证区的滑坡单体划分精度评价结果，本文提出的方法可以将研究区域52.58％的滑坡范围进行正确的单体划分，并且对于面积较大的滑坡体划分的精度较高。
根据本发明构建的滑坡提取方法可以快速有效地提取出滑坡范围，得到较准确的滑坡面积，而本发明提出的滑坡单体划分方法则可以将提取得到的滑坡对象划分成一个个的滑坡单体，避免出现不合理的巨大滑坡区域。联合运用本文构建的滑坡提取与单体划分方法，可对地震或强降雨引起的大范围滑坡灾害进行快速评估，得到滑坡灾害的面积与数量，为灾后应急救援和恢复重建工作提供参考。
附图说明
图1是根据本发明的遥感影像滑坡提取方法的流程图。
图2是根据本发明的遥感影像滑坡单体划分方法的流程图。
图3示出本发明实例中所使用的数据。
图4示出本发明实例的最终滑坡提取结果。
图5示出本发明实例的滑坡连通区域。
图6示出根据本发明实例的滑坡对象融合前后的示意图。
图7示出被滑坡体包围的滑坡对象融合后的示意图。
图8示出滑坡主干与处于边缘的滑坡对象。
图9示出滑坡对象单体划分结果。
图10示出滑坡影像分析对比图。
具体实施方式
下面参照附图对根据本发明的实施例进行详细说明。应当理解，根据本发明的实施例是示意性的而非限定性的。不应将各实施例作为对本发明要求保护的范围的限制。
下面以面向对象的遥感影像滑坡提取方法为例，说明遥感影像滑坡提取方法。本领域技术人员可以理解，其他基于遥感影像的滑坡提取方法例如基于影像上滑坡形态、色调、植被等特征的人工目视解译法以及基于象元光谱信息的图像分类法，均可用于本发明的滑坡单体划分方法。
面向对象分类法不仅可以利用影像的波谱、形态和纹理特征，还可以利用高程数据和其他专题信息，实现滑坡的自动化提取。下面结合图1所示流程图对根据本发明的面向对象的遥感影像滑坡提取方法进行具体说明。
首先，对多光谱影像数据进行多尺度分割，S11。
在研究区域滑坡面积跨度巨大，光谱差异明显的情况下，难于通过一个尺度把滑坡全部分割成单个对象。考虑到滑坡提取时影像分割尺度要尽可能小些，通常选择多个尺度，对影像进行多尺度分割。
随后，对经分割的影像进行滑坡候选对象的提取，S12。
可以通过试错法，选NDVI平均值小于选定阈值的对象为滑坡候选对象。分别合并植被对象和滑坡对象，基本上所有可能的滑坡对象都被归类，与植被对象区分开来。即使有少量处于滑坡边缘的植被象元被划分成了滑坡候选对象，这些象元可以在后续处理中除去。
随后，构造基于光谱、纹理、地形和上下文关系特征的规则，逐步剔除非滑坡对象，获得滑坡对象，S13。
影像中的滑坡由于在光谱特征、几何特征上与滑坡疑似对象特征，例如薄云、河沙、建成区、道路、农田等，容易混淆，在基于对象分类提取滑坡时，往往选择依次剔除非滑坡对象的方法最终提取出真实的滑坡对象。当研究区域范围较大的情况下，不仅地物类型较多，同一类型地物的内部差异也较大。以滑坡为例，其不仅在亮度、湿度等指标上跨度大，而且其形态各异，面积跨度从几百到几百万平方米，坡度跨度从十几到八十度，分布在海拔几百至四千米的范围内。这些都为滑坡的提取增加了难度。本发明根据非滑坡对象的光谱、纹理、地形、上下文关系等构建指标将非滑坡对象提取出来并剔除。各种非滑坡对象的提取指标及提取顺序如表1所示。
对于容易与滑坡对象混淆的非滑坡对象，提取条件设置得越严格，提取出来的非滑坡对象准确性就越高，但与滑坡对象区分度不高的非滑坡对象的遗漏情况就越严重。相反，提取条件设置得越宽松，设置的条件就会越多，虽然遗漏情况得到改善，但非滑坡对象提取的准确性会有所降低。本发明中，采用多次提取法，即每次提取设置不同的条件，或条件的严格程度不同，对某类非滑坡对象进行多次提取，以提取得到准确度高的滑坡对象。
表1非滑坡对象提取特征
对于剔除非滑坡对象后得到的滑坡对象，有的滑坡对象边缘包含着植被象元。随后，采用棋盘分割法，将滑坡对象进行尺度分割。剔除植被，将剩下的滑坡对象融合，作为最终提取滑坡对象，S14。
可以看出，采用如上所述面向对象的滑坡提取方法得到的结果为破碎的滑坡对象。要准确地统计滑坡个数，需要将这些对象划分成一个个的滑坡单体。由于研究区域的滑坡在面积、形状与光谱信息方面都存在巨大的差异，很难运用合适的分割尺度将提取出来的滑坡对象直接分割成滑坡单体。如果直接将这些滑坡对象融合，则会使得许多不同的滑坡体在空间上彼此相接，不利于滑坡个数的统计。
本发明进一步提出了基于面向对象的滑坡单体划分方法。首先通过NDVI标准差阈值使滑坡对象融合形成较大的滑坡对象，然后从目视解译的滑坡结果中总结滑坡面积与滑坡周长的关系、滑坡周长与滑坡共同边界长度的关系，基于此，设置共同边界长度阈值使较大的滑坡对象融合形成滑坡主干。接着，设置共同边界占滑坡周长百分比阈值使得滑坡主干与其周围的滑坡对象融合形成完整的滑坡单体。最后，本发明以目视解译得到的滑坡单体为准，对基于面向对象的滑坡单体划分结果进行验证与精度评价。
本发明认为，彼此相连的不同滑坡，由于空间上紧挨在一起，其滑坡碎屑物相似性较大。研究区域的滑坡绝大多数是震后几个月内发生的新滑坡，由于发生的时间相近，滑坡面有很多相似的地方。这为滑坡的单体划分增加了难度。通过观察，本发明认为滑坡体内部物质间的光谱差异较小，而滑坡体之间的连通区域有较大的光谱差异，突出表现在NDVI上。滑坡的连通区域位于滑坡体的边沿，有限的碎屑物不能把植被完全掩埋，因此，与滑坡主干内部相比，连通区域存在较多NDVI较高的象元。对此，在融合滑坡对象时，通过设置NDVI的标准差阈值，使得滑坡主干与连通区域各自融合成较大的对象。在几何特征上，滑坡连通区域往往较窄，滑坡间的共同边界长度占滑坡体周长的比例较低。利用这一特征，在滑坡单体划分时，通过设置滑坡对象间共同边界的长度阈值、共同边界长度占滑坡周长百分比阈值将不同滑坡体分开。
下面参照如图2示出流程图对根据本发明的面向对象的遥感影像滑坡单体划分方法进行详细说明。
在如图1所示进行基于面向对象的滑坡提取得到最终提取滑坡对象S21后，设置NDVI标准差阈值，将提取的滑坡对象进行多尺度分割融合，S22。对于棋盘分割后剔除了植被象元的对象，面积较小，例如25m2，直接对这些对象进行滑坡单体划分将大大增加计算机的负荷，不利于后续的处理。为了得到最初较为理想尺度下的滑坡对象，对这些小滑坡对象进行多尺度分割融合，分割融合后的滑坡对象仍然非常破碎，需要先划分成较大块的滑坡对象，才能进一步利用共同边界的原则进行单体划分。利用软件中的融合算法(imageobject fusion)，设置融合规划为融合后对象的NDVI标准差不超过设定阈值，使得滑坡主干与连通区域各自融合成较大的对象。由此，滑坡体内部NDVI差异较小的对象能够很好地融合在一起。通过试错法设置NDVI标准差阈值，设置要求是使滑坡体内部对象尽可能融合，而连通区对象不与滑坡体合并。
进一步，对于得到的结果，利用软件中的find enclosed by class算法把被滑坡体包围的滑坡对象融合到一起，以得到更大块的滑坡对象。连通区域由于位于边界，不可能被滑坡对象包围，所以不会被合并。
随后，设置共同边界长度阈值，将共同边界长度大于设定阈值的滑坡对象合并，得到滑坡主干，S23。这个阈值的设定很关键，设得太大，滑坡对象不能融合成一个完整的滑坡，设得太小则会将别的滑坡也包括进来。可以看出，共同边界长度阈值与滑坡的面积和周长有直接的关系。本申请例如采用从滑坡的目视解译结果中提取阈值的确定规则。在最初目视解译出来滑坡对象当中，随机选取满足95％置信度，5％抽样误差的要求的数量的滑坡对象，对于相互连接的不同滑坡体，量测滑坡间最大连通区域的共同边界长度。
在研究区域中滑坡众多，面积从几百平方米到几百万平方米的情况下，大的面积跨度让单体划分时参数的设置变得很困难。如提取大的滑坡，滑坡对象融合时的NDVI标准差阈值要设得小些，同时滑坡对象的共同边界长度阈值要设得大些。而提取小的滑坡时，NDVI标准差阈值要设得大些，共同边界长度阈值要设得小些。因此，本发明根据滑坡面积大小将研究区域的滑坡分为多个等级，逐级进行提取。选择每一级别中间的面积大小，求得其单体划分参数，面积、周长、共同边界、NDVI标准差阈值，作为这一级别滑坡的单体划分参数。例如，对于第一级别滑坡面积，取其中间代表面积，通过幂函数求得其对应的周长，再按共同边界占周长的平均百分比确定其共同边界阈值。
表2示出各个等级滑坡单体划分参数设置。
表2各级滑坡单体划分参数
随后，使用软件中的remove objects算法将共同边界长度大于设定阈值的滑坡对象合并，得到滑坡主干。
此时，滑坡主干周边还可能有许多较小的滑坡对象，其中有一些是属于这个滑坡体的。一类是处于边缘的滑坡对象，因其与滑坡主干有较大的NDVI差异，在第一步的融合时没有包括进来；另一类是对象共同边界长度没有超过设定的阈值，被错误地认为是其他滑坡体。
随后，通过设置滑坡共同边界占周长百分比的阈值将属于滑坡主干的滑坡对象与滑坡主干合并，S24。在上文的统计分析中知道，滑坡体之间的共同边界占滑坡体周长存在最大值。利用这一特点，通过设置滑坡共同边界占周长百分比的阈值可将属于滑坡主干的滑坡对象与滑坡主干合并。考虑到与滑坡主干相连接的滑坡因为NDVI标准差较大等原因，还未融合成一个整体，可适当增大共同边界占滑坡体周长的百分比阈值。利用软件中的remove objects算法，把与滑坡主干的共同边界占其周长超过设定阈值，例如超过设定阈值的5-10％的滑坡对象与滑坡主干合并，由此得到经划分的滑坡单体。
下面以位于汶川地震极重灾区绵竹市的北部及其周边，东经103.94°～104.18°，北纬31.43°～31.66°，南北长约25.7km，东西跨度22.7km，研究区域总面积为583km2为例，具体说明根据本发明的滑坡单体划分方法。研究区域大部分属于绵竹市，其东北角与安县接壤，西北角与茂县相接，而西南角与什邡市毗邻。
本实例使用的是SPOT-5 1A级产品，包括全色与多光谱波段。1A级产品最大程度上保留了卫星原始影像的信息，且提供了星历参数，可用于正向纠正。影像拍摄时间为北京时间2008年10月13日上午11时41分，太阳高度角47.34°。从影像上看，能见度较高，影像右下角存在部分厚云(如图3a)。由太阳高度角不是很大，且研究区域地形崎岖，影像上存在较多阴影。本文利用ENVI 5.1对SPOT-5 1A数据进行辐射校正、几何校正、数据融合等处理得到2.5m的多光谱影像。由于经过了FLAASH大气校正，得到的地表反射率数值范围为0-10000。由于研究区域范围较大，难以获得大比例尺地形图，故采用SRTM空间分辨率为30m的DEM数据。SRTM指的是美国太空署(NASA)与美国国家测绘局(NIMA)共同组织的航天飞机雷达地形测绘任务(Shuttle Radar Topography Mission)。在ENVI(The Environment forVisualizing Images)中运用Topographic Modeling工具通过DEM生成坡度与最大地形曲率数据，导出为图层，参与到滑坡提取过程中(如图3b、图3c)。在ArcGIS中运用水文分析工具对DEM进行洼地填充，将集水面积大于100个象元的汇水区域提取为水系，并利用Strahler算法将提取的水系分为6个等级。通过观察，将第5、6级水系导出为图层(如图3d)，用于滑坡提取时河流对象的剔除。
本实例研究区域滑坡面积跨度巨大，光谱差异明显，考虑到滑坡提取时影像分割尺度要尽可能小，例如选择5个局部峰值中的78为理想的分割尺度，对影像进行多尺度分割。经滑坡候选对象的提取、非滑坡对象的剔除后，采用棋盘分割法，将滑坡对象进行尺度为2的分割。将棋盘分割后平均NDVI≥0.55的对象划分为植被，剩下的则是最终提取的滑坡对象，融合后的结果如图4所示。
本文以目视解译结果为正确的滑坡分布范围，用于验证基于面向对象的滑坡提取结果。在ArcGIS软件中将目视解译结果与滑坡提取结果进行分析。面向对象分类结果与目视解译结果相吻部分为滑坡的正确提取，所以正确识别百分比(Dp)是指正确提取的滑坡面积占目视解译的滑坡面积的百分比。分歧因子(Bf)是指过度提取的滑坡面积占正确提取的滑坡面积的比例。遗漏因子(Mf)是指遗漏提取的滑坡面积占正确提取的滑坡面积的比例。质量百分比(Qp)是指正确提取的滑坡面积占目视解译的滑坡面积与过度提取的滑坡面积之和的百分比。精度评价结果如表3。
表3精度评价
由表3可以看出，根据本发明的面向对象分类的方法能够把研究区域约95％的滑坡提取出来，提取效果较好，分歧因子为0.273。
从总体来看，本实例基于面向对象的滑坡提取的质量百分比为75.30％，提取结果较好。利用构建好的规则集，完成整个研究区域(583km2)的滑坡提取只需约15分钟，说明基于面向对象的方法能够对滑坡进行快速而有效的提取。
随后，在得到的图4所示融合后提取对象的基础上，进行滑坡单体划分。
应当理解，彼此相连的不同滑坡，由于空间上紧挨在一起，其滑坡碎屑物相似性较大。而且研究区域的滑坡绝大多数是震后几个月内发生的新滑坡，由于发生的时间相近，滑坡面有很多相似的地方。这为滑坡的单体划分增加了难度。通过观察，可以看出滑坡体内部物质间的光谱差异较小，而滑坡体之间的连通区域有较大的光谱差异，突出表现在NDVI上。如图5中的椭圆部分，滑坡的连通区域位于滑坡体的边沿，有限的碎屑物不能把植被完全掩埋，因此，与滑坡主干内部相比，连通区域存在较多NDVI较高的象元。对此，本实例在融合滑坡对象时，通过设置NDVI的标准差阈值，使得滑坡主干与连通区域各自融合成较大的对象。在几何特征上，滑坡连通区域往往较窄(如图5椭圆部分)，滑坡间的共同边界长度占滑坡体周长的比例较低。利用这一特征，在滑坡单体划分时，通过设置滑坡对象间共同边界的长度阈值、共同边界长度占滑坡周长百分比阈值将不同滑坡体分开。
对于棋盘分割后剔除了植被象元的对象，面积较小(25m2)，直接对这些对象进行滑坡单体划分将大大增加计算机的负荷，不利于后续的处理。为了得到最初较为理想尺度下的滑坡对象，对这些小滑坡对象进行多尺度分割融合，参数的设置与最初的影像分割相同(分割尺度＝78，形状系数＝0.1，颜色系数＝0.9，紧致度＝0.5，光滑度＝0.5)，结果如图6a。分割融合后的滑坡对象仍然非常破碎，需要先划分成较大块的滑坡对象，才能进一步利用共同边界的原则进行单体划分。利用软件中的融合算法(image object fusion)，设置融合规划为融合后对象的NDVI标准差不超过某个阈值，使得滑坡主干与连通区域各自融合成较大的对象。如图6b，滑坡体内部NDVI差异较小的对象能够很好地融合在一起。通过试错法设置NDVI标准差阈值，设置要求是使滑坡体内部对象尽可能融合，而连通区对象不与滑坡体合并。
对于得到的结果，本文利用软件中的find enclosed by class算法把被滑坡体包围的滑坡对象融合到一起，以得到更大块的滑坡对象。连通区域由于位于边界，不可能被滑坡对象包围，所以不会被合并，结果如图7。
需要设置共同边界长度阈值把大块的滑坡对象合并成单个的滑坡体。很明显，共同边界长度阈值与滑坡的面积和周长有直接的关系。本实例从滑坡的目视解译结果中提取阈值的确定规则。在最初目视解译出来的2185个滑坡当中随机选取351个滑坡(满足95％置信度，5％抽样误差的要求)，对于相互连接的不同滑坡体，量测滑坡间最大连通区域的共同边界长度。在抽取的351个滑坡体当中，有81个滑坡体出现与周边滑坡体相连的现象，因此测得81个共同边界的长度。在ArcGIS中使用Calculate Geometry工具计算抽取的滑坡体面积与周长。分析这些滑坡体的面积与其周长的关系，发现滑坡体面积与其周长的关系可用幂函数y＝0.505x0.5117(R2＝0.9382)来刻画。统计81个滑坡体的共同边界与其周长的关系，绝大多数(87％)滑坡体间的共同边界小于滑坡体周长的20％，而最大值为29.86％，最小值为1.23％，平均值为10.75％。
本实例中，研究区域滑坡众多，面积从几百平方米到几百万平方米。如此大的面积跨度让单体划分时参数的设置变得很困难。如提取大的滑坡，滑坡对象融合时的NDVI标准差阈值要设得小些，同时滑坡对象的共同边界长度阈值要设得大些。而提取小的滑坡时，NDVI标准差阈值要设得大些，共同边界长度阈值要设得小些。本实例中，根据滑坡面积大小将研究区域的滑坡分为七个等级，逐级进行提取。选择每一级别中间的面积大小，求得其单体划分参数，作为这一级别滑坡的单体划分参数。如第一级别滑坡面积为1百万至1千万平方米，则其中间代表面积为5百万平方米，通过幂函数求得其对应的周长为5.505*(5000000)0.5117＝14744m，再按共同边界占周长的平均百分比确定其共同边界阈值14744*10.75％＝1585m。各个等级滑坡单体划分参数设置如表2。使用软件中的remove objects算法将共同边界长度大于其阈值的滑坡对象合并，得到滑坡主干，如图8中的中部大部分区域。
根据图8，滑坡主干周边还有许多较小的滑坡对象，其中有一些是属于这个滑坡体的。一类是处于边缘的滑坡对象，因其与滑坡主干有较大的NDVI差异，在第一步的融合时没有包括进来(如图8中的中上部椭圆)；另一类是对象共同边界长度没有超过设定的阈值，被错误地认为是其他滑坡体(如图8中的圆圈)。
根据上文的统计分析可以看出，滑坡体之间的共同边界占滑坡体周长的最大值为29.86％，也就是说不超过滑坡体周长的30％。利用这一特点，通过设置滑坡共同边界占周长百分比的阈值将属于滑坡主干的滑坡对象与主滑坡合并。考虑到与滑坡主干相连接的滑坡因为NDVI标准差较大等原因，还未融合成一个整体(如图9中的中部和中下部椭圆)，所以适当增大共同边界占滑坡体周长的百分比阈值。利用软件中的remove objects算法，把与滑坡主干的共同边界占其周长超过35％的滑坡对象与滑坡主干合并，结果如图9。可以看出，滑坡单体的边缘相对光滑，形状指数较小。若某个滑坡体形状指数特别大，说明其边缘非常破碎，它很可能是由几个滑坡连在一起形成的。通过试验，设置形状指数阈值为4.5，小于这个阈值的为滑坡单体，大于这一阈值的则参与下一级别的滑坡单体划分。下一级别滑坡的提取重复以上步骤，只是在设置NDVI标准差阈值时要适当逐级增大，因为这可以使无法参与上一级滑坡融合的对象参与到这一级的滑坡融合当中。本文中NDVI标准差阈值随着滑坡面积的减小按一定的步长增大。
到划分第七级面积小于1000m2的小型滑坡时，有少数属于前面级别，面积较大的滑坡因为通不过形状指数检验而留在这里。对此，在进行单体划分的前两步(NDVI标准差阈值融合、被包围的滑坡对象合并)后，只利用前面滑坡等级的共同边界长度阈值将这些对象合并，按面积划分到各个滑坡等级中，不再考虑别的因素。对于最后面积小于1000m2的滑坡对象则直接融合，划分为第七级滑坡。
结果验证与精度评价
图10示出流域单体划分与其他结果对比结果。图10a为影像图，图10b为融合结果，图10c为提取的滑坡，图10d为划分的单体图。整个研究区域滑坡对象直接融合得到的结果如图10b，而滑坡单体划分的最终结果如图10d，两者存在巨大的差异。从图上可以看出，直接融合可能造成整个流域的滑坡体连成一片，形成巨大的滑坡区域。而单体划分则可以有效避免这种情况的发生。对比单体划分的结果与直接融合的结果(如表4)，单体划分后面积大于1000m2的滑坡个数约为直接融合的2倍，而滑坡的平均面积则是其二分之一。直接融合的滑坡面积标准差超过23万平方米，而单体划分后的标准差只有约8万平方米。这些数据说明滑坡的单体划分能有效地将大量空间相连的滑坡区域进行划分，避免出现不合理的巨大滑坡。
表4滑坡单体划分结果与直接融合结果对比
本发明提出的基于面向对象的滑坡单体划分方法，可以快速有效地对得到的滑坡对象进行单体划分，对于准确统计滑坡数量有重要意义。利用本文提出的方法构建好规则集，可以在3分钟内完成研究区域112.46km2滑坡对象的单体划分，得到10090个滑坡单体。与滑坡对象直接融合的结果相比，滑坡单体划分得到的结果可以避免形成不合理的巨大滑坡区域。
根据三个验证区的滑坡单体划分精度评价结果，本文提出的方法可以对研究区域52.58％的滑坡范围进行正确的单体划分，滑坡单体划分后得到的滑坡个数约为直接融合的2倍，平均滑坡面积约为直接融合的一半，而滑坡面积的标准差则只有直接融合的36％。
本发明在接近现实情况的较大范围研究区域内，详细探讨了基于面向对象的遥感影像滑坡提取方法，并提出了滑坡单体划分方法。结果表明，本发明构建的滑坡提取方法可以快速有效地提取出滑坡范围，得到较准确的滑坡面积，而本发明提出的滑坡单体划分方法则可以有效地将滑坡对象划分成一个个的滑坡单体，避免出现不合理的大滑坡区域。联合运用本发明提出的滑坡提取与单体划分方法，可对地震或强降雨引起的大范围滑坡灾害进行快速评估，得到滑坡灾害的面积与数量，为灾后应急救援和恢复重建工作提供参考。
显而易见的是，在不脱离本发明的精神和范围的情况下，在阅读前述公开内容之后，本发明的各种其它变型和调整对于本领域技术人员将是显而易见的，并且意图是所有这样的变型和调整在所附权利要求的范围内。
Claims (10)
1.一种遥感影像滑坡单体划分方法，其特征在于，该方法包括：
以研究区域的多光谱影像数据和DEM数据为基础，进行滑坡提取，得到滑坡对象；
设置NDVI标准差阈值，利用融合算法将所得到的滑坡对象融合；
利用滑坡间的共同边界长度阈值使经融合的滑坡对象进一步融合得到滑坡主干；以及
利用共同边界占滑坡周长百分比阈值使得滑坡主干与其周围的滑坡对象融合，得到滑坡单体。
2.如权利要求1所述的遥感影像滑坡单体划分方法，其特征在于，该方法进一步包括根据滑坡面积大小将研究区域的滑坡对象分为多个等级，分别对各等级的滑坡对象设置NDVI标准差阈值、共同边界长度阈值和共同边界占滑坡周长百分比阈值，进行融合。
3.如权利要求1所述的遥感影像滑坡单体划分方法，其特征在于，所述设置NDVI标准差阈值包括：
通过试错法设置NDVI标准差阈值，设置要求是使滑坡体滑坡对象尽可能融合，而连通区内滑坡对象不与滑坡体合并。
4.如权利要求1所述的遥感影像滑坡单体划分方法，其特征在于，按滑坡对象共同边界占周长的平均百分比确定其共同边界长度阈值。
5.如权利要求1所述的遥感影像滑坡单体划分方法，其特征在于，滑坡对象共同边界占周长的百分比阈值大于该滑坡对象共同边界占周长的平均百分比。
6.如权利要求1所述的遥感影像滑坡单体划分方法，其特征在于，所述以研究区域的多光谱影像数据和DEM数据为基础，进行滑坡提取得到滑坡对象的步骤进一步包括：
对多光谱影像数据进行多尺度分割；
对经分割的影像进行滑坡候选对象的提取，
构造基于光谱、纹理、地形和上下文关系特征的规则，逐步剔除非滑坡对象，获得滑坡对象。
7.如权利要求6所述的遥感影像滑坡单体划分方法，其特征在于，所述对经分割的影像进行滑坡候选对象的提取的步骤包括利用NDVI滑坡提取阈值区分植被对象和滑坡对象，提取滑坡候选对象。
8.如权利要求6所述的遥感影像滑坡单体划分方法，其特征在于，所述逐步剔除非滑坡对象的步骤进一步包括，按照积雪、云、阴影、河、湖泊、人造物的顺序剔除非滑坡对象。
9.如权利要求8所述的遥感影像滑坡单体划分方法，其特征在于，该方法进一步包括对得到的滑坡对象进行棋盘分割，以剔除小的植被对象，得到提取的滑坡对象。
10.如权利要求1所述的遥感影像滑坡单体划分方法，其特征在于，该方法进一步包括，利用Find Enclosed by Class算法，将被滑坡体包围的滑坡对象融合在一起。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710393741.4A CN107255516B (zh) | 2017-05-27 | 2017-05-27 | 一种遥感影像滑坡单体划分方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710393741.4A CN107255516B (zh) | 2017-05-27 | 2017-05-27 | 一种遥感影像滑坡单体划分方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107255516A CN107255516A (zh) | 2017-10-17 |
| CN107255516B true CN107255516B (zh) | 2018-11-06 |
Family
ID=60027512
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710393741.4A Active CN107255516B (zh) | 2017-05-27 | 2017-05-27 | 一种遥感影像滑坡单体划分方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107255516B (zh) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108827880B (zh) * | 2018-04-23 | 2021-04-16 | 吉林大学 | 基于多光谱影像和ndvi时间序列的地表覆盖变化检测方法 |
| CN108875615B (zh) * | 2018-06-07 | 2021-04-30 | 中国石油天然气股份有限公司 | 沉积区域遥感识别方法、装置、电子设备及存储介质 |
| CN111626269B (zh) * | 2020-07-07 | 2021-08-27 | 中国科学院空天信息创新研究院 | 一种实用的大空间范围滑坡提取方法 |
| CN111882573B (zh) * | 2020-07-31 | 2023-08-18 | 北京师范大学 | 一种基于高分辨率影像数据的耕地地块提取方法及系统 |
| CN112419495B (zh) * | 2020-10-26 | 2022-11-15 | 天津大学 | 基于多尺度dem空间模型的高程点自动提取方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2947727B2 (ja) * | 1995-04-18 | 1999-09-13 | 明星電気株式会社 | 地下湧水の検知方法及び地下湧水監視装置 |
| CN105989322B (zh) * | 2015-01-27 | 2019-07-05 | 同济大学 | 一种基于高分辨率遥感影像的多指标融合滑坡检测方法 |
-
2017
- 2017-05-27 CN CN201710393741.4A patent/CN107255516B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN107255516A (zh) | 2017-10-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107255516B (zh) | 一种遥感影像滑坡单体划分方法 | |
| Pradhan | Groundwater potential zonation for basaltic watersheds using satellite remote sensing data and GIS techniques | |
| Saha et al. | Land cover classification using IRS LISS III image and DEM in a rugged terrain: a case study in Himalayas | |
| Ismail et al. | Satellite data classification accuracy assessment based from reference dataset | |
| Schneevoigt et al. | Detecting Alpine landforms from remotely sensed imagery. A pilot study in the Bavarian Alps | |
| TW200929067A (en) | 3D image detecting, editing and rebuilding system | |
| Kasprzak et al. | UAV and SfM in detailed geomorphological mapping of granite tors: an example of Starościńskie Skały (Sudetes, SW Poland) | |
| Mohamed et al. | Impacts of soil sealing on potential agriculture in Egypt using remote sensing and GIS techniques | |
| Bar Massada et al. | Automated segmentation of vegetation structure units in a Mediterranean landscape | |
| CN105447274A (zh) | 一种利用面向对象分类技术对中等分辨率遥感图像进行滨海湿地制图的方法 | |
| Clark et al. | Landscape analysis using multi-scale segmentation and object-oriented classification | |
| CN106485718A (zh) | 一种过火迹地识别方法及装置 | |
| Geerling et al. | Mapping river floodplain ecotopes by segmentation of spectral (CASI) and structural (LiDAR) remote sensing data | |
| Broussard III et al. | Quantifying vegetation and landscape metrics with hyperspatial unmanned aircraft system imagery in a coastal oligohaline marsh | |
| Franklin | Land cover stratification using Landsat Thematic Mapper data in Sahelian and Sudanian woodland and wooded grassland | |
| Singh et al. | Effect of sensor modelling methods on computation of 3-D coordinates from Cartosat-1 stereo data | |
| Jeziorska et al. | Overland flow analysis using time series of sUAS-derived elevation models | |
| Quattrochi et al. | Fractal characterization of multitemporal remote sensing data | |
| Härmä et al. | The production of finnish CORINE land cover 2000 classification | |
| Yimer et al. | Seasonal effect on the accuracy of Land use/Land cover classification in the Bilate Sub-basin, Abaya-Chamo Basin, Rift valley Lakes Basin of Ethiopia. | |
| CN109885808A (zh) | 一种近地表气象要素计算方法 | |
| Tian et al. | A process-oriented method for rapid acquisition of canopy height model from RGB point cloud in semiarid region | |
| Göksel et al. | Land Use and Land Cover Changes Using Spot 5 Pansharpen Images; A Case Study in Akdeniz District, Mersin-Turkey | |
| Saba et al. | Co-seismic landslides automatic detection on regional scale with sub-pixel analysis of multi temporal high resolution optical images: Application to southwest of Port Au Prince, Haiti | |
| Nguyen-Van-Anh et al. | The influence of satellite image spatial resolution on mapping land use/land cover: A case study of Ho Chi Minh City, Vietnam |
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 |