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

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

【技术实现步骤摘要】
基于遥感影像的泥石流灾害识别及频率计算方法


[0001]本专利技术涉及基于摄影测量的地质灾害调查领域，具体涉及基于遥感影像的泥石流灾害识别及频率计算方法。

技术介绍

[0002]泥石流是我国山区常见的地质灾害类型之一，对山区经济发展和人民安全构成严重威胁。由于山区泥石流沟数量众多、分布广泛，且每条泥石流沟发生泥石流灾害的频次均有所差异，精准掌握泥石流灾害频率信息，对于科学评估和防治泥石流灾害具有重要意义。
[0003]当前，泥石流灾害的发生频率信息主要依靠两种方式获取：一、现场调查；二、文献调查。
[0004]其中，现场调查，由于泥石流分布范围广阔，且部分区域存在通行障碍，人员现场调查费时费力，仅能够解决一部分位于公路沿线临近的泥石流沟的灾害频率信息获取问题，此外，相关信息只是被参与现场调查的相关人员掌握，难于共享，不利于区域联防联控。
[0005]文献调查，则是从相关报道、新闻等文献资料中进行查询，能够快速获得一些泥石流沟的灾害发生频次信息，但这些信息能够覆盖的泥石流沟数量有限，仅限于少量影响重大的灾害事件，且频次信息受限于截止至资料发表时间之前，后续的灾害频次信息难以持续更新，对广域范围内泥石流灾害评估和防治工作的信息贡献较为乏力。
[0006]因此，目前，尚未出现针对广域范围内泥石流灾害频率信息计算的深入研究。

技术实现思路

[0007]本专利技术所要解决的技术问题是：提出一种基于遥感影像的泥石流灾害识别及频率计算方法，解决传统方式获取泥石流灾害频率信息存在的效率低、覆盖面不足、频次信息难以更新的问题。
[0008]本专利技术解决上述技术问题采用的技术方案是：基于遥感影像的泥石流灾害识别方法，包括以下步骤：S1、数据准备：根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集；根据目标泥石流流域所在区域的数字高程模型，提取目标泥石流流域的边界范围，并基于流域内坡度的变化将流域划分为子区域；S2、根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数，然后，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数；S3、基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像
素数据集；S4、基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件。
[0009]具体的，步骤S1中，根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集，包括：A1、从遥感影像数据集U
L
中，按自然年度分别提取待识别年度及其前后两年的遥感影像，并构建各年度的年度数据集U
LY
，下标Y表示年度；A2、逐景判定年度数据集U
LY
中各影像数据I
ij
的所属月份，所述I
ij
表示第j年的第i景的遥感影像数据；然后，对连续三年的遥感影像数据，根据夏半年及冬半年的划分进行重构，构建待识别年度的夏半年的前后两个冬半年的数据集U
DY
，下标Y表示年度。
[0010]进一步的，步骤A2中，在根据夏半年及冬半年的划分进行重构后，首先，将待录入冬半年数据集U
DY
的各影像数据I
ij
，逐个计算其云层覆盖率，然后，仅将满足云层覆盖率条件的影像数据录入数据集，完成两个冬半年数据集U
DY
的构建。
[0011]具体的，提取的遥感影像数据，其光谱波段包括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
的光谱信息，逐像素的进行云层像素判定，并将满足如下组合条件的像素认定为云层像素：
其中，表示第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
的云层像素的数量，按如下公式计算其云层覆盖率：其中，C
ij
为第j年第i景影像数据的云层覆盖率，N
ij
‑
Cloud
为影像数据I
ij
中云层像素的数量，N
ij
‑
All
为影像数据I
ij
中的总像素数量。
[0012]具体的，步骤S1中，根据目标泥石流流域所在区域的数字高程模型，采用ArcGIS软件，提取目标泥石流流域的边界范围。
[0013]具体的，步骤S1中，基于流域内坡度的变化将流域划分为子区域，包括：B1、根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集；B2、遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点；B3、基于分区突变点，计算子区域的分界线，并将边界范围内的目标泥石流流域划分为子区域；B4、通过匹配数字高程模型和卫星影像数据的位置信息，提取卫星影像数据中与目标泥石流流域的边界范围、分界线及主河道对应像素的坐标。
[0014]最优的，步骤B1中，根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集，包括：B11、遍历边界范围内目标泥石流流域每个点位的高程值，并将最高高程点e
max
作为目标泥石流流域的主沟道起始点D1；B12、将主沟道起始点D1的空间位置及高程值录入主沟道点集，并将主沟道起始点D1作为中心点D
C
的初始输入；B13、根据输入的中心点D
C
，遍历其周围的八个点位，将其中高程值最低的点位识别为主沟道点D
k
，将D
k
的空间位置及高程值录入主沟道点集，下标k为主沟道点的序号；B14、判定是否到达目标泥石流流域的边界，若是，则完成流域的主沟道的识别；否则，将主沟道点D
k
作为新的中心点D
C
，返回步骤B13。
[0015]最优的，目标泥石流流域的子区域包括形成区、流通区和堆积区，步骤B2中，遍本文档来自技高网...

【技术保护点】

【技术特征摘要】
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
的光谱信息，逐像素的进行云层像素判定，并将满足如下组合条件的像素认定为云层像素：
其中，表示第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
的云层像素的数量，按如下公式计算其云层覆盖率：其中，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
作为目标泥石流流域的主沟道起始点D1；B12、将主沟道起始点D1的空间位置及高程值录入主沟道点集，并将主沟道起始点D1作为中心点D
C
的初始输入；B13、根据输入的中心点D
C
，遍历其周围的八个点位，将其中高程值最低的点位识别为主沟道点D
k
，将D
k
的空间位置及高程值录入主沟道点集，下标k为主沟道点的序号；B14、判定是否到达目标泥石流流域的边界，若是，则完成流域的主沟道的识别；否则，将主沟道点D
k
作为新的中心点D
C
，返回步骤B13。8.如权利要求6所述的基于遥感影像的泥石流灾害识别方法，其特征在于，目标泥石流流域的子区域包括形成区、流通区和堆积区，步骤B2中，遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点，包括：B21、从主沟起始点D1开始，依序计算主沟道点集中相邻两点间的纵坡比降：其中，S
k
为主沟道点D
k
与D
k+1
两点间的纵坡比降，
∆
h
k
为主沟道点D
k
与D
k+1
两点间的高程差值，
∆
l
k
为主沟道点D
k
与D
k+1
两点间的水平距离，下标k为主沟道点的序号；B22、将主沟道进行分段，并基于各分段所包含主沟道点的纵坡比降，计算分段的平均纵坡比降：其中，为第q个分段的平均纵坡比降，S
qs
为第q个分段所包含的第s个主沟道点的纵坡比降，r为分段所包含主沟道点的数量；B23、从主沟起始点D1开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于或等于50%时，将连续多个分段中的首个分段作为形成区和流通区的划界分段；然后，从形成区和流通区的划界分段开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于1%时，将连续多个分段中的首个分段作为流通区和堆积区的划界分段；B24、基于形成区和流通区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点D
min
‑
A1
；基于流通区和堆积区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点D
min
‑
A2
，以D
min
‑
A1
和D
min
‑
A2
作为分区突变点。9.如权利要求8所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤B24中，从形成区和流通区的划界分段中，提取分段所包含主沟道点中纵坡比降最
小的主沟道点D<...

【专利技术属性】
技术研发人员：张建强，贾洋，张莉莉，葛永刚，黄宇，明仔洋，
申请(专利权)人：四川省公路规划勘察设计研究院有限公司，
类型：发明
国别省市：