CN115272874B - 基于遥感影像的泥石流灾害识别及频率计算方法 - Google Patents

基于遥感影像的泥石流灾害识别及频率计算方法 Download PDF

Info

Publication number
CN115272874B
CN115272874B CN202211187713.4A CN202211187713A CN115272874B CN 115272874 B CN115272874 B CN 115272874B CN 202211187713 A CN202211187713 A CN 202211187713A CN 115272874 B CN115272874 B CN 115272874B
Authority
CN
China
Prior art keywords
debris flow
year
main channel
pixel
remote sensing
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
CN202211187713.4A
Other languages
English (en)
Other versions
CN115272874A (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.)
Institute of Mountain Hazards and Environment IMHE of CAS
Sichuan Highway Planning Survey and Design Institute Ltd
Original Assignee
Institute of Mountain Hazards and Environment IMHE of CAS
Sichuan Highway Planning Survey and Design Institute Ltd
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 Institute of Mountain Hazards and Environment IMHE of CAS, Sichuan Highway Planning Survey and Design Institute Ltd filed Critical Institute of Mountain Hazards and Environment IMHE of CAS
Priority to CN202211187713.4A priority Critical patent/CN115272874B/zh
Publication of CN115272874A publication Critical patent/CN115272874A/zh
Application granted granted Critical
Publication of CN115272874B publication Critical patent/CN115272874B/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
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/13—Satellite images
    • 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/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/77—Processing image or video features in feature spaces; using data integration or data reduction, e.g. principal component analysis [PCA] or independent component analysis [ICA] or self-organising maps [SOM]; Blind source separation
    • G06V10/774—Generating sets of training patterns; Bootstrap methods, e.g. bagging or boosting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes
    • G06V20/188—Vegetation

Landscapes

  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Multimedia (AREA)
  • General Health & Medical Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Computing Systems (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Evolutionary Computation (AREA)
  • Databases & Information Systems (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Astronomy & Astrophysics (AREA)
  • Remote Sensing (AREA)
  • Image Processing (AREA)

Abstract

本发明涉及基于摄影测量的地质灾害调查领域，提供了一种基于遥感影像的泥石流灾害识别及频率计算方法，解决传统方式获取泥石流灾害频率信息存在的效率低、覆盖面不足、频次信息难以更新的问题。本发明的识别方法，首选筛选待识别年度及其前后两年的遥感影像，同时基于高程数据提取泥石流流域边界和分区信息；然后，根据各子区域在待识别年度前后两个冬半年的影像数据中相应波段的TOA值，计算像素级的归一化植被指数和物源扰动变化检测指数，进而识别各子区域物源的变化情况，最终，基于变化并结合判定条件，识别是否发生灾害事件；而频率计算方法，首先，基于识别方法识别统计周期内的灾害发生次数，然后计算泥石流的发生频率。

Description

基于遥感影像的泥石流灾害识别及频率计算方法
技术领域
本发明涉及基于摄影测量的地质灾害调查领域，具体涉及基于遥感影像的泥石流灾害识别及频率计算方法。
背景技术
泥石流是我国山区常见的地质灾害类型之一，对山区经济发展和人民安全构成严重威胁。由于山区泥石流沟数量众多、分布广泛，且每条泥石流沟发生泥石流灾害的频次均有所差异，精准掌握泥石流灾害频率信息，对于科学评估和防治泥石流灾害具有重要意义。
当前，泥石流灾害的发生频率信息主要依靠两种方式获取：一、现场调查；二、文献调查。
其中，现场调查，由于泥石流分布范围广阔，且部分区域存在通行障碍，人员现场调查费时费力，仅能够解决一部分位于公路沿线临近的泥石流沟的灾害频率信息获取问题，此外，相关信息只是被参与现场调查的相关人员掌握，难于共享，不利于区域联防联控。
文献调查，则是从相关报道、新闻等文献资料中进行查询，能够快速获得一些泥石流沟的灾害发生频次信息，但这些信息能够覆盖的泥石流沟数量有限，仅限于少量影响重大的灾害事件，且频次信息受限于截止至资料发表时间之前，后续的灾害频次信息难以持续更新，对广域范围内泥石流灾害评估和防治工作的信息贡献较为乏力。
因此，目前，尚未出现针对广域范围内泥石流灾害频率信息计算的深入研究。
发明内容
本发明所要解决的技术问题是：提出一种基于遥感影像的泥石流灾害识别及频率计算方法，解决传统方式获取泥石流灾害频率信息存在的效率低、覆盖面不足、频次信息难以更新的问题。
本发明解决上述技术问题采用的技术方案是：
基于遥感影像的泥石流灾害识别方法，包括以下步骤：
S1、数据准备：
根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集；
根据目标泥石流流域所在区域的数字高程模型，提取目标泥石流流域的边界范围，并基于流域内坡度的变化将流域划分为子区域；
S2、根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数，然后，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数；
S3、基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集；
S4、基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件。
具体的，步骤S1中，根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集，包括：
A1、从遥感影像数据集U L 中，按自然年度分别提取待识别年度及其前后两年的遥感影像，并构建各年度的年度数据集U LY ，下标Y表示年度；
A2、逐景判定年度数据集U LY 中各影像数据I ij 的所属月份，所述I ij 表示第j年的第i景的遥感影像数据；然后，对连续三年的遥感影像数据，根据夏半年及冬半年的划分进行重构，构建待识别年度的夏半年的前后两个冬半年的数据集U DY ，下标Y表示年度。
进一步的，步骤A2中，在根据夏半年及冬半年的划分进行重构后，首先，将待录入冬半年数据集U DY 的各影像数据I ij ，逐个计算其云层覆盖率，然后，仅将满足云层覆盖率条件的影像数据录入数据集，完成两个冬半年数据集U DY 的构建。
具体的，提取的遥感影像数据，其光谱波段包括0.45~0.52μm、0.52~0.60μm、0.63~0.69μm、0.76~0.90μm、1.55~1.75μm、10.40~12.50μm和2.08~2.35μm七个波段；
在步骤A2中，仅将满足云层覆盖率<1%的影像数据录入数据集，并按如下步骤计算云层覆盖率：
A21、根据影像数据I ij 的光谱信息，逐像素的进行云层像素判定，并将满足如下组合条件的像素认定为云层像素：
Figure DEST_PATH_IMAGE001
其中，
Figure DEST_PATH_IMAGE002
表示第j年第i景的影像数据的第m个像素的光谱信息中，在第n个波段的TOA值；所述B1~B7依次对应0.45~0.52μm、0.52~0.60μm、0.63~0.69μm、0.76~0.90μm、1.55~1.75μm、10.40~12.50μm和2.08~2.35μm七个波段；
A22、根据影像数据I ij 的云层像素的数量，按如下公式计算其云层覆盖率：
Figure DEST_PATH_IMAGE003
其中，C ij 为第j年第i景影像数据的云层覆盖率，N ij-Cloud 为影像数据I ij 中云层像素的数量，N ij-All 为影像数据I ij 中的总像素数量。
具体的，步骤S1中，根据目标泥石流流域所在区域的数字高程模型，采用ArcGIS软件，提取目标泥石流流域的边界范围。
具体的，步骤S1中，基于流域内坡度的变化将流域划分为子区域，包括：
B1、根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集；
B2、遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点；
B3、基于分区突变点，计算子区域的分界线，并将边界范围内的目标泥石流流域划分为子区域；
B4、通过匹配数字高程模型和卫星影像数据的位置信息，提取卫星影像数据中与目标泥石流流域的边界范围、分界线及主河道对应像素的坐标。
最优的，步骤B1中，根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集，包括：
B11、遍历边界范围内目标泥石流流域每个点位的高程值，并将最高高程点e max 作为目标泥石流流域的主沟道起始点D 1 ；
B12、将主沟道起始点D 1 的空间位置及高程值录入主沟道点集，并将主沟道起始点D 1 作为中心点D C 的初始输入；
B13、根据输入的中心点D C ，遍历其周围的八个点位，将其中高程值最低的点位识别为主沟道点D k ，将D k 的空间位置及高程值录入主沟道点集，下标k为主沟道点的序号；
B14、判定是否到达目标泥石流流域的边界，若是，则完成流域的主沟道的识别；否则，将主沟道点D k 作为新的中心点D C ，返回步骤B13。
最优的，目标泥石流流域的子区域包括形成区、流通区和堆积区，步骤B2中，遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点，包括：
B21、从主沟起始点D 1 开始，依序计算主沟道点集中相邻两点间的纵坡比降：
Figure DEST_PATH_IMAGE004
其中，S k 为主沟道点D k 与D k+1 两点间的纵坡比降，∆h k 为主沟道点D k 与D k+1 两点间的高程差值，∆l k 为主沟道点D k 与D k+1 两点间的水平距离，下标k为主沟道点的序号；
B22、将主沟道进行分段，并基于各分段所包含主沟道点的纵坡比降，计算分段的平均纵坡比降：
Figure DEST_PATH_IMAGE005
其中，
Figure DEST_PATH_IMAGE006
为第q个分段的平均纵坡比降，S qs 为第q个分段所包含的第s个主沟道点的纵坡比降，r为分段所包含主沟道点的数量；
B23、从主沟起始点D 1 开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于或等于50%时，将连续多个分段中的首个分段作为形成区和流通区的划界分段；然后，从形成区和流通区的划界分段开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于1%时，将连续多个分段中的首个分段作为流通区和堆积区的划界分段；
B24、基于形成区和流通区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A1 ；基于流通区和堆积区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A2 ，以D min-A1 和D min-A2 作为分区突变点。
最优的，步骤B24中，从形成区和流通区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A1 和纵坡比降最大的主沟道点D max-A1 ；从流通区和堆积区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A2 和纵坡比降最大的主沟道点D max-A2 ；且提取的主沟道点D min-A1 、D max-A1 、D min-A2 和D max-A2 的排序应满足：
Figure DEST_PATH_IMAGE007
其中，D 1 为主沟道的起始点，D K 为主沟道的终点；
以D min-A1 和D min-A2 作为分区突变点，并将D max-A1 和D max-A2 作为辅助计算点；
在步骤B3中，在数字高程模型的水平投影内，将通过点D min-A1 且垂直于D min-A1 和D max-A1 连线的直线作为形成区和流通区的分界线；将通过点D min-A2 且垂直于D min-A2 和D max-A2 连线的直线作为流通区和堆积区的分界线。
最优的，步骤B22中，将主沟道按长度80~120米进行分段，步骤B23中，所述连续多个分段包括：连续三个分段。
具体的，步骤S2，首先，按下式，根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数：
Figure DEST_PATH_IMAGE008
其中，
Figure DEST_PATH_IMAGE009
和
Figure DEST_PATH_IMAGE010
为第j个冬半年的第i景影像数据的第m个像素在0.63~0.69μm波段和0.76~0.90μm波段的TOA值；
然后，按下式，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数：
Figure DEST_PATH_IMAGE011
其中，j表示在前的冬半年， j+1表示在后的冬半年。
具体的，步骤S3中，基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集，包括：
S31、基于步骤S2中，由前后两个冬半年的第一景影像数据，计算获得的归一化植被指数NDVI m1j+1 和物源扰动变化检测指数Index m1j+1 ，按如下条件逐像素进行判定，并将满足条件的像素录入新增物源待定像素数据集P PD ：
Figure DEST_PATH_IMAGE012
其中，Thred NDVI 为归一化植被指数的判定阈值，Thred Index 为物源扰动变化检测指数的判定阈值；
S32、基于步骤S2中，由前后两个冬半年的第一景以外的其他景影像数据，计算获得的归一化植被指数NDVI mkj+1 和物源扰动变化检测指数Index mkj+1 ，分别将各景对应数据中与数据集P PD 各像素相同坐标的像素的数据，按如下条件逐景逐像素进行判定：
Figure DEST_PATH_IMAGE013
其中，k∈(1,M],M=min{N j ,N j+1 }，N j 为第j个冬半年影像数据的景数；
根据在各景中均满足条件的像素的坐标，对数据集P PD 中的像素进行筛选，形成新增物源像素数据集；
S33、根据新增物源像素数据集所包含各像素的坐标，按目标泥石流流域各子区域间的分界线进性重构，获得各子区域的新增物源像素数据集。
具体的，目标泥石流流域的子区域包括形成区、流通区和堆积区；
步骤S4中，基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件，包括：
S41、将流通区新增物源像素数据集P LT 中各像素与主河道的最小距离，分别与预设的流通区空间距离阈值Thred LT 进行比较，若小于阈值Thred LT ，则判定该像素为流通区泥石流活动有效堆积物源像素；然后，基于判定结果，构建流通区的有效堆积物源像素数据集P LT-final ；
将堆积区新增物源像素数据集P DJ 中各像素与主河道的最小距离，分别与预设的堆积区空间距离阈值Thred DJ 进行比较，若存在小于阈值Thred DJ 的像素，则判定新增物源像素数据集P DJ 所包含像素均为堆积区泥石流活动有效堆积物源像素，否则，判定新增物源像素数据集P DJ 所包含像素均为堆积区泥石流活动无效堆积物源像素；然后，基于判定结果，构建堆积区的有效堆积物源像素数据集P DJ-final ；
S42、若流通区的有效堆积物源像素数据集P LT-final 的像素数量大于或等于阈值Thred LT-pixel ，或堆积区的有效堆积物源像素数据集P DJ-final 的像素数量大于或等于阈值Thred DJ-pixel ，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
进一步的，在步骤S42完成后，还包括：
S43、若形成区新增物源像素数据集P XC 的像素数量大于或等于阈值Thred XC-pixel ，且流通区的有效堆积物源像素数据集P LT-final 或堆积区的有效堆积物源像素数据集P DJ-final 的像素数量大于或等于阈值Thred LT-DJ ，则，进一步统计待识别年度连续15日最大累计降雨量P Acc15 ，若待识别年度的P Acc15 超过待识别年度前3年同期的平均水平，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
另外，本发明还提供了一种基于遥感影像的泥石流灾害频率计算方法，包括：
首先，根据如上所述的基于遥感影像的泥石流灾害识别方法，对指定的统计周期内的年度，逐年识别目标泥石流流域是否发生泥石流；
然后，按如下公式计算目标泥石流流域的灾害发生频率F debrisFlow ：
Figure DEST_PATH_IMAGE014
其中，N event 为统计周期内目标泥石流流域发生泥石流灾害事件的次数，N year 为统计周期的年份数。
本发明的有益效果是：
本发明基于遥感影像，利用像元光谱差异化特征计算，能够实现灾害发生前泥石流沟道物源区松散物源自动识别，及灾害发生后在泥石流沟道及沟口堆积物质的自动提取；结合地貌形态特征因素，通过时间序列对比计算，精准提取时间周期内堆积物变化情况，实现灾害发生、规模体量的量化判定，进而实现泥石流沟灾害发生频率的快速计算。该方法适用于广域范围的大量泥石流沟的灾害识别及频率信息提取，革新了传统的泥石流频率信息获取方式。
附图说明
图1为本发明基于遥感影像的泥石流灾害识别方法的实现框架图。
具体实施方式
本发明旨在提出一种适用于广域范围内的基于遥感影像的泥石流灾害识别方法，及基于其的泥石流灾害频率计算方法，以解决传统方式获取泥石流灾害频率信息存在的效率低、覆盖面不足、频次信息难以更新的问题。
本发明的基于遥感影像的泥石流灾害识别方法，包括：
S1、数据准备：
根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集；
根据目标泥石流流域所在区域的数字高程模型，提取目标泥石流流域的边界范围，并基于流域内坡度的变化将流域划分为子区域；
S2、根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数，然后，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数；
S3、基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集；
S4、基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件。
下面结合实施例，对本发明的方案作进一步的描述。
实施例：
如图1所示，本实施例的广域范围泥石流灾害的识别及频率计算方法，包括以下步骤：
S1、数据准备，其包括地形数据和遥感影像数据的准备。
A、遥感影像数据的准备：
根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集。
由于泥石流主要由降雨导致，通常发生在夏半年，因此，泥石流灾害的发生，可通过其前后两个冬半年的影像对比加以识别。在北半球，夏半年，通常指春分日到秋分日之间的这段时间；冬半年，通常指秋季10月经冬季到来年春季3月的这段时间，为了方便统一计算，在本实施例中，冬半年具体指当年11月、12月和来年的1月和2月这段时间。
具体的，包括：
A1、从遥感影像数据集U L 中，按自然年度分别提取待识别年度及其前后两年的遥感影像，并构建各年度的年度数据集U LY ，下标Y表示年度；
A2、逐景判定年度数据集U LY 中各影像数据I ij 的所属月份，所述I ij 表示第j年的第i景的遥感影像数据；然后，对连续三年的遥感影像数据，根据夏半年及冬半年的划分进行重构，构建待识别年度的夏半年的前后两个冬半年的数据集U DY ，下标Y表示年度。
如待识别年度为2003年，则提取2002、2003和2004连续三年的数据，并构建2002年11月~2003年2月以及2003年11月~2004年2月两个冬半年的数据集。上述的景是卫星遥感影像数据的单位，1景即卫星拍摄的1张照片，同一区域卫星在某一天过境只拍摄1张照片，即1景影像，时间区间内，不同重访周期对应相应区域的影像数据景数=时间跨度总天数/重访周期天数。
在高原及山区，冬半年通常属旱季，通常无云，但为了排除云层的干扰，进一步的，步骤A2中，在根据夏半年及冬半年的划分进行重构后，首先，将待录入冬半年数据集U DY 的各影像数据I ij ，逐个计算其云层覆盖率，然后，仅将满足云层覆盖率条件的影像数据录入数据集，完成两个冬半年数据集U DY 的构建。云层覆盖率的计算可以采用现有任意云检测算法，其选择与遥感数据所包含的光谱信息有关。
在本实施例中，提取的遥感影像数据，其光谱波段包括：波长0.45~0.52μm的蓝绿光光谱、波长0.52~0.60μm的绿光光谱、波长0.63~0.69μm的红光光谱、波长0.76~0.90μm的近红外光光谱、波长1.55~1.75μm的中红外光光谱、波长10.40~12.50μm的热红外光光谱和波长2.08~2.35μm的短波红外光光谱，七个波段。因此，实施例中，采用Fmask云检测算法。同时，为了严格排除云层影响，实施例中，仅将满足云层覆盖率<1%的影像数据录入数据集。
具体的，在步骤A2中，并按如下步骤计算云层覆盖率：
A21、根据影像数据I ij 的光谱信息，逐像素的进行云层像素判定，并将满足如下组合条件的像素认定为云层像素：
Figure 221783DEST_PATH_IMAGE001
其中，
Figure 935661DEST_PATH_IMAGE002
表示第j年第i景的影像数据的第m个像素的光谱信息中，在第n个波段的TOA值；所述B1~B7依次对应0.45~0.52μm、0.52~0.60μm、0.63~0.69μm、0.76~0.90μm、1.55~1.75μm、10.40~12.50μm和2.08~2.35μm七个波段；
A22、根据影像数据I ij 的云层像素的数量，按如下公式计算其云层覆盖率：
Figure 512136DEST_PATH_IMAGE003
其中，C ij 为第j年第i景影像数据的云层覆盖率，N ij-Cloud 为影像数据I ij 中云层像素的数量，N ij-All 为影像数据I ij 中的总像素数量。
B、地形数据的准备：
首先，根据目标泥石流流域所在区域的数字高程模型，提取目标泥石流流域的边界范围，并基于流域内坡度的变化将流域划分为子区域。目标泥石流流域的边界范围可以采用现有任意的算法，在本实施例中，采用ArcGIS软件，提取目标泥石流流域的边界范围。
然后，基于流域内坡度的变化将流域划分为子区域，具体包括：
B1、根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集。
流域中，一定是由高往低流动，因此，可以采用逐点比较的方式来识别主沟道，但其同时受到分辨率及微地形的影响，在分辨率较高的情况下，也可以采用多点平均高程或坡度的方式来克服微地形的影响，
在本实施例中，数字高程模型的分辨率为30米，因此，采用逐点比较的方式，具体的，步骤B1包括：
B11、遍历边界范围内目标泥石流流域每个点位的高程值，并将最高高程点e max 作为目标泥石流流域的主沟道起始点D 1 ；
B12、将主沟道起始点D 1 的空间位置及高程值录入主沟道点集，并将主沟道起始点D 1 作为中心点D C 的初始输入；
B13、根据输入的中心点D C ，遍历其周围的八个点位，将其中高程值最低的点位识别为主沟道点D k ，将D k 的空间位置及高程值录入主沟道点集，下标k为主沟道点的序号；
B14、判定是否到达目标泥石流流域的边界，若是，则完成流域的主沟道的识别；否则，将主沟道点D k 作为新的中心点D C ，返回步骤B13。
B2、遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点。
子区域的划分可以根据需要进行划分，而要实现灾害发生前泥石流沟道物源区松散物源自动识别，及灾害发生后在泥石流沟道及沟口堆积物质自动提取，因此，最好的，将目标泥石流流域的子区域包括形成区、流通区和堆积区，具体的，步骤B2包括：
B21、从主沟起始点D 1 开始，依序计算主沟道点集中相邻两点间的纵坡比降：
Figure 468197DEST_PATH_IMAGE004
其中，S k 为主沟道点D k 与D k+1 两点间的纵坡比降，∆h k 为主沟道点D k 与D k+1 两点间的高程差值，∆l k 为主沟道点D k 与D k+1 两点间的水平距离，下标k为主沟道点的序号；
B22、将主沟道进行分段，并基于各分段所包含主沟道点的纵坡比降，计算分段的平均纵坡比降：
Figure 813728DEST_PATH_IMAGE005
其中，
Figure 432928DEST_PATH_IMAGE006
为第q个分段的平均纵坡比降，S qs 为第q个分段所包含的第s个主沟道点的纵坡比降，r为分段所包含主沟道点的数量；
B23、从主沟起始点D 1 开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于或等于50%时，将连续多个分段中的首个分段作为形成区和流通区的划界分段；然后，从形成区和流通区的划界分段开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于1%时，将连续多个分段中的首个分段作为流通区和堆积区的划界分段；
B24、基于形成区和流通区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A1 ；基于流通区和堆积区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A2 ，以D min-A1 和D min-A2 作为分区突变点。
上述的判定条件：平均纵坡比降均小于或等于50%、平均纵坡比降均小于1%，是根据川西地区已有的泥石流灾害样本统计获得的形成区、流通区和堆积区之间的分界地形判定阈值条件。当然，针对其他的区域，可以根据该区域的具体情况，如已有的泥石流灾害样本，进行设置；同时，若对流域子区域的划分采用其他方式，阈值条件应与之相匹配。
上述的，将主沟道进行分段，同时基于连续多个分段判定，是为了避免微观地形，比如巨石或者断崖等的影响，而分段的长度以及连续多个分段的数量则与数字高程模型的分辨率有关。最优的，步骤B22中，将主沟道按长度80~120米进行分段，步骤B23中，所述连续多个分段包括：连续三个分段。在实施例中，数字高程模型的分辨率为30米，每段主沟道包括3个点位，由于识别的主沟道点连线呈折线状，因此，3个点位，能覆盖约100米的范围。
B3、基于分区突变点，计算子区域的分界线，并将边界范围内的目标泥石流流域划分为子区域。分界线也即经过分区突变点的线，并基于该分界线将流域划分为子区域，其可根据需要进行设定，比如，按主沟道两侧坡地的坡度设置斜线进行分界。
在本实施例中，为了方便计算，具体的，包括：
步骤B24中，从形成区和流通区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A1 和纵坡比降最大的主沟道点D max-A1 ；从流通区和堆积区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A2 和纵坡比降最大的主沟道点D max-A2 ；且提取的主沟道点D min-A1 、D max-A1 、D min-A2 和D max-A2 的排序应满足：
Figure 496699DEST_PATH_IMAGE007
其中，D 1 为主沟道的起始点，D K 为主沟道的终点；上述的(D 1 ，D min-A1 ］表示D max-A1 位于主沟起始点D 1 和D min-A1 之间，其中，圆括号表示不包括该点，方括号表示包括该点，其他区间同理。
以D min-A1 和D min-A2 作为分区突变点，并将D max-A1 和D max-A2 作为辅助计算点；
在步骤B3中，在数字高程模型的水平投影内，将通过点D min-A1 且垂直于D min-A1 和D max-A1 连线的直线作为形成区和流通区的分界线；将通过点D min-A2 且垂直于D min-A2 和D max-A2 连线的直线作为流通区和堆积区的分界线。
B4、通过匹配数字高程模型和卫星影像数据的位置信息，提取卫星影像数据中与目标泥石流流域的边界范围、分界线及主河道对应像素的坐标。
S2、首先，根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数。
归一化植被指数，英文缩写为NDVI，能反映出植物冠层的背景影响，如土壤、潮湿地面、雪、枯叶、粗糙度等信息，且与植被覆盖有关。因此，基于前后两个冬半年的归一化植被指数，能获得各像素的变化情况，也即后述的像素级的物源扰动变化检测指数。
归一化植被指数的计算，如下式所述：
Figure DEST_PATH_IMAGE015
其中，
Figure DEST_PATH_IMAGE016
和
Figure DEST_PATH_IMAGE017
为第j个冬半年的第i景影像数据的第m个像素在0.63~0.69μm波段和0.76~0.90μm波段的TOA值；
然后，按下式，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数：
Figure 118435DEST_PATH_IMAGE011
其中，j表示在前的冬半年， j+1表示在后的冬半年。
S3、基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集。
为了降低计算量，在本实施例中，具体的，步骤S3包括：
S31、基于步骤S2中，由前后两个冬半年的第一景影像数据，计算获得的归一化植被指数NDVI m1j+1 和物源扰动变化检测指数Index m1j+1 ，按如下条件逐像素进行判定，并将满足条件的像素录入新增物源待定像素数据集P PD ：
Figure 584052DEST_PATH_IMAGE012
其中，Thred NDVI 为归一化植被指数的判定阈值，Thred Index 为物源扰动变化检测指数的判定阈值；
S32、基于步骤S2中，由前后两个冬半年的第一景以外的其他景影像数据，计算获得的归一化植被指数NDVI mkj+1 和物源扰动变化检测指数Index mkj+1 ，分别将各景对应数据中与数据集P PD 各像素相同坐标的像素的数据，按如下条件逐景逐像素进行判定：
Figure 374153DEST_PATH_IMAGE013
其中，k∈(1,M],M=min{N j ,N j+1 }，N j 为第j个冬半年影像数据的景数；
根据在各景中均满足条件的像素的坐标，对数据集P PD 中的像素进行筛选，形成新增物源像素数据集；
S33、根据新增物源像素数据集所包含各像素的坐标，按目标泥石流流域各子区域间的分界线进性重构，获得各子区域的新增物源像素数据集。
在本实施例中，上述的阈值Thred NDVI 和Thred Index ，根据川西地区已有灾害样本统计分析得到，Thred NDVI 设定为0.05，Thred Index 设定为0.45。当然，针对其他的区域，可以根据该区域的具体情况，如已有的泥石流灾害样本，进行设置。
S4、基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件。
灾害阈值判定条件，主要根据已有灾害样本统计分析得到，并同子区域的划分有关。而要实现灾害发生前泥石流沟道物源区松散物源自动识别，也即对形成区物源变化的识别；及灾害发生后在泥石流沟道及沟口堆积物质自动提取，也即，流通区和堆积区物源变化的识别，因此，最好的，将目标泥石流流域的子区域包括形成区、流通区和堆积区。
因此，具体的，步骤S4包括：
S41、将流通区新增物源像素数据集P LT 中各像素与主河道的最小距离，分别与预设的流通区空间距离阈值Thred LT 进行比较，若小于阈值Thred LT ，则判定该像素为流通区泥石流活动有效堆积物源像素；然后，基于判定结果，构建流通区的有效堆积物源像素数据集P LT-final ；
将堆积区新增物源像素数据集P DJ 中各像素与主河道的最小距离，分别与预设的堆积区空间距离阈值Thred DJ 进行比较，若存在小于阈值Thred DJ 的像素，则判定新增物源像素数据集P DJ 所包含像素均为堆积区泥石流活动有效堆积物源像素，否则，判定新增物源像素数据集P DJ 所包含像素均为堆积区泥石流活动无效堆积物源像素；然后，基于判定结果，构建堆积区的有效堆积物源像素数据集P DJ-final ；
S42、若流通区的有效堆积物源像素数据集P LT-final 的像素数量大于或等于阈值Thred LT-pixel ，或堆积区的有效堆积物源像素数据集P DJ-final 的像素数量大于或等于阈值Thred DJ-pixel ，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
S43、若形成区新增物源像素数据集P XC 的像素数量大于或等于阈值Thred XC-pixel ，且流通区的有效堆积物源像素数据集P LT-final 或堆积区的有效堆积物源像素数据集P DJ-final 的像素数量大于或等于阈值Thred LT-DJ ，则，进一步统计待识别年度连续15日最大累计降雨量P Acc15 ，若待识别年度的P Acc15 超过待识别年度前3年同期的平均水平，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
同样的，根据对川西地区已有灾害样本的统计分析，Thred LT-pixel 设定为5，Thred DJ-pixel 设定为10； Thred XC-pixel 设定为5，Thred LT-DJ 设定为1。Thred LT 设定为60米，Thred DJ 设定为30米。而上述的各像素与主河道的最小距离，也即该像素与主河道对应像素间距离的最小值。
连续15日最大累计降雨量P Acc15 ，也即一年中，以15日为单位，计算获得的累计降雨量的最大值，具体计算见下式：
Figure DEST_PATH_IMAGE018
上式中，下标mj表示第j年的m日。
进一步的，基于上述的识别方法，能够实现基于遥感影像的泥石流灾害频率计算方法，包括：
首先，根据前述的基于遥感影像的泥石流灾害识别方法，对指定的统计周期内的年度，逐年识别目标泥石流流域是否发生泥石流；
然后，按如下公式计算目标泥石流流域的灾害发生频率F debrisFlow ：
Figure 226353DEST_PATH_IMAGE014
其中，N event 为统计周期内目标泥石流流域发生泥石流灾害事件的次数，N year 为统计周期的年份数。
尽管这里参照本发明的实施例对本发明进行了描述，上述实施例仅为本发明较佳的实施方式，本发明的实施方式并不受上述实施例的限制，应该理解，本领域技术人员可以设计出很多其他的修改和实施方式，这些修改和实施方式将落在本申请公开的原则范围和精神之内。

Claims (15)

1.基于遥感影像的泥石流灾害识别方法，其特征在于，包括以下步骤：
S1、数据准备：
根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集；
根据目标泥石流流域所在区域的数字高程模型，提取目标泥石流流域的边界范围，并基于流域内坡度的变化将流域划分为子区域；
S2、根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数，然后，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数；
S3、基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集；
S4、基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件。
2.如权利要求1所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤S1中，根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集，包括：
A1、从遥感影像数据集U L 中，按自然年度分别提取待识别年度及其前后两年的遥感影像，并构建各年度的年度数据集U LY ，下标Y表示年度；
A2、逐景判定年度数据集U LY 中各影像数据I ij 的所属月份，所述I ij 表示第j年的第i景的遥感影像数据；然后，对连续三年的遥感影像数据，根据夏半年及冬半年的划分进行重构，构建待识别年度的夏半年的前后两个冬半年的数据集U DY ，下标Y表示年度。
3.如权利要求2所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤A2中，在根据夏半年及冬半年的划分进行重构后，首先，将待录入冬半年数据集U DY 的各影像数据I ij ，逐个计算其云层覆盖率，然后，仅将满足云层覆盖率条件的影像数据录入数据集，完成两个冬半年数据集U DY 的构建。
4.如权利要求3所述的基于遥感影像的泥石流灾害识别方法，其特征在于，提取的遥感影像数据，其光谱波段包括0.45~0.52μm、0.52~0.60μm、0.63~0.69μm、0.76~0.90μm、1.55~1.75μm、10.40~12.50μm和2.08~2.35μm七个波段；
在步骤A2中，仅将满足云层覆盖率<1%的影像数据录入数据集，并按如下步骤计算云层覆盖率：
A21、根据影像数据I ij 的光谱信息，逐像素的进行云层像素判定，并将满足如下组合条件的像素认定为云层像素：
Figure 96768DEST_PATH_IMAGE001
其中，
Figure 279487DEST_PATH_IMAGE002
表示第j年第i景的影像数据的第m个像素的光谱信息中，在第n个波段的TOA值；B1~B7依次对应0.45~0.52μm、0.52~0.60μm、0.63~0.69μm、0.76~0.90μm、1.55~1.75μm、10.40~12.50μm和2.08~2.35μm七个波段；
A22、根据影像数据I ij 的云层像素的数量，按如下公式计算其云层覆盖率：
Figure 387121DEST_PATH_IMAGE003
其中，C ij 为第j年第i景影像数据的云层覆盖率，N ij-Cloud 为影像数据I ij 中云层像素的数量，N ij-All 为影像数据I ij 中的总像素数量。
5.如权利要求1所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤S1中，根据目标泥石流流域所在区域的数字高程模型，采用ArcGIS软件，提取目标泥石流流域的边界范围。
6.如权利要求1所述的基于遥感影像的泥石流灾害识别方法，其特征在于，
步骤S1中，基于流域内坡度的变化将流域划分为子区域，包括：
B1、根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集；
B2、遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点；
B3、基于分区突变点，计算子区域的分界线，并将边界范围内的目标泥石流流域划分为子区域；
B4、通过匹配数字高程模型和卫星影像数据的位置信息，提取卫星影像数据中与目标泥石流流域的边界范围、分界线及主河道对应像素的坐标。
7.如权利要求6所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤B1中，根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集，包括：
B11、遍历边界范围内目标泥石流流域每个点位的高程值，并将最高高程点e max 作为目标泥石流流域的主沟道起始点D 1 ；
B12、将主沟道起始点D 1 的空间位置及高程值录入主沟道点集，并将主沟道起始点D 1 作为中心点D C 的初始输入；
B13、根据输入的中心点D C ，遍历其周围的八个点位，将其中高程值最低的点位识别为主沟道点D k ，将D k 的空间位置及高程值录入主沟道点集，下标k为主沟道点的序号；
B14、判定是否到达目标泥石流流域的边界，若是，则完成流域的主沟道的识别；否则，将主沟道点D k 作为新的中心点D C ，返回步骤B13。
8.如权利要求6所述的基于遥感影像的泥石流灾害识别方法，其特征在于，目标泥石流流域的子区域包括形成区、流通区和堆积区，步骤B2中，遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点，包括：
B21、从主沟起始点D 1 开始，依序计算主沟道点集中相邻两点间的纵坡比降：
Figure 844647DEST_PATH_IMAGE004
其中，S k 为主沟道点D k 与D k+1 两点间的纵坡比降，∆h k 为主沟道点D k 与D k+1 两点间的高程差值，∆l k 为主沟道点D k 与D k+1 两点间的水平距离，下标k为主沟道点的序号；
B22、将主沟道进行分段，并基于各分段所包含主沟道点的纵坡比降，计算分段的平均纵坡比降：
Figure 924598DEST_PATH_IMAGE005
其中，
Figure 310842DEST_PATH_IMAGE006
为第q个分段的平均纵坡比降，S qs 为第q个分段所包含的第s个主沟道点的纵坡比降，r为分段所包含主沟道点的数量；
B23、从主沟起始点D 1 开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于或等于50%时，将连续多个分段中的首个分段作为形成区和流通区的划界分段；然后，从形成区和流通区的划界分段开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于1%时，将连续多个分段中的首个分段作为流通区和堆积区的划界分段；
B24、基于形成区和流通区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A1 ；基于流通区和堆积区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A2 ，以D min-A1 和D min-A2 作为分区突变点。
9.如权利要求8所述的基于遥感影像的泥石流灾害识别方法，其特征在于，
步骤B24中，从形成区和流通区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A1 和纵坡比降最大的主沟道点D max-A1 ；从流通区和堆积区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点D min-A2 和纵坡比降最大的主沟道点D max-A2 ；且提取的主沟道点D min-A1 、D max-A1 、D min-A2 和D max-A2 的排序应满足：
Figure 905772DEST_PATH_IMAGE007
其中，D 1 为主沟道的起始点，D K 为主沟道的终点；
以D min-A1 和D min-A2 作为分区突变点，并将D max-A1 和D max-A2 作为辅助计算点；
在步骤B3中，在数字高程模型的水平投影内，将通过点D min-A1 且垂直于D min-A1 和D max-A1 连线的直线作为形成区和流通区的分界线；将通过点D min-A2 且垂直于D min-A2 和D max-A2 连线的直线作为流通区和堆积区的分界线。
10.如权利要求8或9任一项所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤B22中，将主沟道按长度80~120米进行分段，步骤B23中，所述连续多个分段包括：连续三个分段。
11.如权利要求1~9任一项所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤S2，首先，按下式，根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数：
Figure 370251DEST_PATH_IMAGE008
其中，
Figure 367026DEST_PATH_IMAGE009
和
Figure 891548DEST_PATH_IMAGE010
为第j个冬半年的第i景影像数据的第m个像素在0.63~0.69μm波段和0.76~0.90μm波段的TOA值；
然后，按下式，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数：
Figure 415851DEST_PATH_IMAGE011
其中，j表示在前的冬半年， j+1表示在后的冬半年。
12.如权利要求11所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤S3中，基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集，包括：
S31、基于步骤S2中，由前后两个冬半年的第一景影像数据，计算获得的归一化植被指数NDVI m1j+1 和物源扰动变化检测指数Index m1j+1 ，按如下条件逐像素进行判定，并将满足条件的像素录入新增物源待定像素数据集P PD ：
Figure 480759DEST_PATH_IMAGE012
其中，Thred NDVI 为归一化植被指数的判定阈值，Thred Index 为物源扰动变化检测指数的判定阈值；
S32、基于步骤S2中，由前后两个冬半年的第一景以外的其他景影像数据，计算获得的归一化植被指数NDVI mkj+1 和物源扰动变化检测指数Index mkj+1 ，分别将各景对应数据中与数据集P PD 各像素相同坐标的像素的数据，按如下条件逐景逐像素进行判定：
Figure 66462DEST_PATH_IMAGE013
其中，k∈(1,M],M=min{N j ,N j+1 }，N j 为第j个冬半年影像数据的景数；
根据在各景中均满足条件的像素的坐标，对数据集P PD 中的像素进行筛选，形成新增物源像素数据集；
S33、根据新增物源像素数据集所包含各像素的坐标，按目标泥石流流域各子区域间的分界线进性重构，获得各子区域的新增物源像素数据集。
13.如权利要求1~9任一项所述的基于遥感影像的泥石流灾害识别方法，其特征在于，目标泥石流流域的子区域包括形成区、流通区和堆积区；
步骤S4中，基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件，包括：
S41、将流通区新增物源像素数据集P LT 中各像素与主河道的最小距离，分别与预设的流通区空间距离阈值Thred LT 进行比较，若小于阈值Thred LT ，则判定该像素为流通区泥石流活动有效堆积物源像素；然后，基于判定结果，构建流通区的有效堆积物源像素数据集P LT-final ；
将堆积区新增物源像素数据集P DJ 中各像素与主河道的最小距离，分别与预设的堆积区空间距离阈值Thred DJ 进行比较，若存在小于阈值Thred DJ 的像素，则判定新增物源像素数据集P DJ 所包含像素均为堆积区泥石流活动有效堆积物源像素，否则，判定新增物源像素数据集P DJ 所包含像素均为堆积区泥石流活动无效堆积物源像素；然后，基于判定结果，构建堆积区的有效堆积物源像素数据集P DJ-final ；
S42、若流通区的有效堆积物源像素数据集P LT-final 的像素数量大于或等于阈值Thred LT-pixel ，或堆积区的有效堆积物源像素数据集P DJ-final 的像素数量大于或等于阈值Thred DJ-pixel ，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
14.如权利要求13所述的基于遥感影像的泥石流灾害识别方法，其特征在于，在步骤S42完成后，还包括：
S43、若形成区新增物源像素数据集P XC 的像素数量大于或等于阈值Thred XC-pixel ，且流通区的有效堆积物源像素数据集P LT-final 或堆积区的有效堆积物源像素数据集P DJ-final 的像素数量大于或等于阈值Thred LT-DJ ，则，进一步统计待识别年度连续15日最大累计降雨量P Acc15 ，若待识别年度的P Acc15 超过待识别年度前3年同期的平均水平，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
15.基于遥感影像的泥石流灾害频率计算方法，其特征在于，
首先，根据权利要求1~14任一项所述的基于遥感影像的泥石流灾害识别方法，对指定的统计周期内的年度，逐年识别目标泥石流流域是否发生泥石流；
然后，按如下公式计算目标泥石流流域的灾害发生频率F debrisFlow ：
Figure 558623DEST_PATH_IMAGE014
其中，N event 为统计周期内目标泥石流流域发生泥石流灾害事件的次数，N year 为统计周期的年份数。
CN202211187713.4A 2022-09-28 2022-09-28 基于遥感影像的泥石流灾害识别及频率计算方法 Active CN115272874B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211187713.4A CN115272874B (zh) 2022-09-28 2022-09-28 基于遥感影像的泥石流灾害识别及频率计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211187713.4A CN115272874B (zh) 2022-09-28 2022-09-28 基于遥感影像的泥石流灾害识别及频率计算方法

Publications (2)

Publication Number Publication Date
CN115272874A CN115272874A (zh) 2022-11-01
CN115272874B true CN115272874B (zh) 2022-12-09

Family

ID=83756371

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211187713.4A Active CN115272874B (zh) 2022-09-28 2022-09-28 基于遥感影像的泥石流灾害识别及频率计算方法

Country Status (1)

Country Link
CN (1) CN115272874B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116091881B (zh) * 2023-02-14 2023-07-11 安徽星太宇科技有限公司 一种基于多源数据融合的遥感信息管理系统
CN115908954B (zh) * 2023-03-01 2023-07-28 四川省公路规划勘察设计研究院有限公司 基于人工智能的地质灾害隐患识别系统、方法及电子设备
CN116127787B (zh) * 2023-04-07 2023-06-23 中国科学院、水利部成都山地灾害与环境研究所 火烧烈度-高程积分方法及火后泥石流易发性评估方法
CN117809260B (zh) * 2024-02-23 2024-05-14 中国地质调查局水文地质环境地质调查中心 泥石流物源启动的识别方法、装置以及电子设备

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004244947A (ja) * 2003-02-14 2004-09-02 Sabo Jisuberi Gijutsu Center 土砂災害の氾濫解析システム
CN103530516A (zh) * 2013-10-14 2014-01-22 成都理工大学 一种强震区泥石流隐患点快速识别方法
CN109767409A (zh) * 2018-11-28 2019-05-17 中国科学院遥感与数字地球研究所 基于遥感影像的滑坡变化检测方法、存储介质和电子设备
CN110765650A (zh) * 2019-11-12 2020-02-07 中国科学院、水利部成都山地灾害与环境研究所 一种泥石流体积含沙量的测算方法
CN111832430A (zh) * 2020-06-23 2020-10-27 长江水利委员会长江科学院 山洪灾害高光谱遥感影像识别方法及其识别系统
CN113538861A (zh) * 2021-09-16 2021-10-22 山东金科星机电股份有限公司 基于矿产地质勘查的地质灾害信息管理系统

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004244947A (ja) * 2003-02-14 2004-09-02 Sabo Jisuberi Gijutsu Center 土砂災害の氾濫解析システム
CN103530516A (zh) * 2013-10-14 2014-01-22 成都理工大学 一种强震区泥石流隐患点快速识别方法
CN109767409A (zh) * 2018-11-28 2019-05-17 中国科学院遥感与数字地球研究所 基于遥感影像的滑坡变化检测方法、存储介质和电子设备
CN110765650A (zh) * 2019-11-12 2020-02-07 中国科学院、水利部成都山地灾害与环境研究所 一种泥石流体积含沙量的测算方法
CN111832430A (zh) * 2020-06-23 2020-10-27 长江水利委员会长江科学院 山洪灾害高光谱遥感影像识别方法及其识别系统
CN113538861A (zh) * 2021-09-16 2021-10-22 山东金科星机电股份有限公司 基于矿产地质勘查的地质灾害信息管理系统

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
基于时间序列遥感影像的滑坡检测方法;虎振兴 等;《航天返回与遥感》;20180415;第39卷(第2期);104-114 *
空间高分辨率数字正射影像黄土区地质灾害解译特征与应用: 以甘肃天水麦积区幅为例;田尤 等;《现代地质》;20161031;第30卷(第5期);1161-1169 *
西藏洛隆县巴曲冰湖溃决型泥石流演进过程模拟研究;刘波 等;《水文地质工程地质》;20190930;第48卷(第5期);150-160 *
高寒高海拔山原区沟谷型泥石流成因与特征———以四川省雅江县祝桑景区为例;倪化勇 等;《水土保持通报》;20130228;第33卷(第1期);211-215 *

Also Published As

Publication number Publication date
CN115272874A (zh) 2022-11-01

Similar Documents

Publication Publication Date Title
CN115272874B (zh) 基于遥感影像的泥石流灾害识别及频率计算方法
Brandolini et al. Response of terraced slopes to a very intense rainfall event and relationships with land abandonment: A case study from Cinque Terre (Italy)
Zhang et al. Urban built-up land change detection with road density and spectral information from multi-temporal Landsat TM data
Van den Berg et al. Detection, quantification and monitoring of Prosopis in the Northern Cape Province of South Africa using remote sensing and GIS
CN109635731B (zh) 一种识别有效耕地的方法及装置、存储介质及处理器
Dufour et al. Image utilisation for the study and management of riparian vegetation: overview and applications
Chouari Wetland land cover change detection using multitemporal Landsat data: A case study of the Al-Asfar wetland, Kingdom of Saudi Arabia
Dashpurev et al. A cost-effective method to monitor vegetation changes in steppes ecosystems: A case study on remote sensing of fire and infrastructure effects in eastern Mongolia
CN114881454A (zh) 山体滑坡与崩塌灾害分类方法与电子设备
Ochieng et al. Use of remote sensing and geographical information system (GIS) for salinity assessment of Vaal-Harts irrigation scheme, South Africa
Gardiner et al. A quantitative appraisal of woody shrub encroachment in western New South Wales.
CN114881457B (zh) 基于决策树的山体滑坡与崩塌灾害分类方法与电子设备
Kinga The spatio-temporal analysis of impervious surfaces in Cluj-Napoca, Romania
Saha et al. Investigation of land-use change and groundwater–surface water interaction in the Kiskatinaw River watershed, northeastern British Columbia (parts of NTS 093P/01,/02,/07–/10)
Pani et al. Delineation and monitoring of gullied and ravinous lands in a part of lower Chambal Valley, India, using remote sensing and GIS
Han et al. Water distribution based on SAR and optical data to improve hazard mapping
Shih Satellite data and geographic information system for land use classification
Azzam Land suitability evaluation for cultivation of some soils in western desert of egypt, el-minya governorate using gis and remote sensing
Owe et al. Improved classification of small-scale urban watersheds using thematic mapper simulator data
Sader et al. Time-series tropical forest change detection: a visual and quantitative approach
Andrade An exploratory study on heads up photo interpretation of aerial photography as a method for mapping drainage tiles
Baučić et al. Open Geospatial Data for Urban Green Areas
Zia et al. Impacts of Urbanization on Green Spaces of the Densely Populated City of Karachi, Pakistan-An Analysis of 8 Years of Data for Estimating Land Cover Changes
Rawat et al. Aeolian sand affected wasteland monitoring using multi-temporal remotely sensed imagery: A case study of Sirsa district of Haryana, India
CN115049931A (zh) 一种基于多光谱成像仪的水土保持措施智能化提取和评价方法

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