技术特征：
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
min-a1
和纵坡比降最大的主沟道点d
max-a1
；从流通区和堆积区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点d
min-a2
和纵坡比降最大的主沟道点d
max-a2
；且提取的主沟道点d
min-a1
、d
max-a1
、d
min-a2
和d
max-a2
的排序应满足：其中，d1为主沟道的起始点，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值，获取各遥感影像的像素级的归一化植被指数：其中，和 为第j个冬半年的第i景影像数据的第m个像素在0.63~0.69μm波段和0.76~0.90μm波段的toa值；然后，按下式，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数：其中，j表示在前的冬半年， j+1表示在后的冬半年。12.如权利要求11所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤s3中，基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集，包括：s31、基于步骤s2中，由前后两个冬半年的第一景影像数据，计算获得的归一化植被指数ndvi
m1j+1
和物源扰动变化检测指数index
m1j+1
，按如下条件逐像素进行判定，并将满足条件的像素录入新增物源待定像素数据集p
pd
：
其中，thred
ndvi
为归一化植被指数的判定阈值，thred
index
为物源扰动变化检测指数的判定阈值；s32、基于步骤s2中，由前后两个冬半年的第一景以外的其他景影像数据，计算获得的归一化植被指数ndvi
mkj+1
和物源扰动变化检测指数index
mkj+1
，分别将各景对应数据中与数据集p
pd
各像素相同坐标的像素的数据，按如下条件逐景逐像素进行判定：其中，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
：其中，n
event
为统计周期内目标泥石流流域发生泥石流灾害事件的次数，n
year
为统计周期的年份数。
技术总结
本发明涉及基于摄影测量的地质灾害调查领域，提供了一种基于遥感影像的泥石流灾害识别及频率计算方法，解决传统方式获取泥石流灾害频率信息存在的效率低、覆盖面不足、频次信息难以更新的问题。本发明的识别方法，首选筛选待识别年度及其前后两年的遥感影像，同时基于高程数据提取泥石流流域边界和分区信息；然后，根据各子区域在待识别年度前后两个冬半年的影像数据中相应波段的TOA值，计算像素级的归一化植被指数和物源扰动变化检测指数，进而识别各子区域物源的变化情况，最终，基于变化并结合判定条件，识别是否发生灾害事件；而频率计算方法，首先，基于识别方法识别统计周期内的灾害发生次数，然后计算泥石流的发生频率。率。率。
技术研发人员：张建强 贾洋 张莉莉 葛永刚 黄宇 明仔洋
受保护的技术使用者：四川省公路规划勘察设计研究院有限公司
技术研发日：2022.09.28
技术公布日：2022/11/2