本发明涉及基于摄影测量的地质灾害调查领域，提供了一种基于遥感影像的泥石流灾害识别及频率计算方法，解决传统方式获取泥石流灾害频率信息存在的效率低、覆盖面不足、频次信息难以更新的问题。本发明的识别方法，首选筛选待识别年度及其前后两年的遥感影像，同时基于高程数据提取泥石流流域边界和分区信息；然后，根据各子区域在待识别年度前后两个冬半年的影像数据中相应波段的TOA值，计算像素级的归一化植被指数和物源扰动变化检测指数，进而识别各子区域物源的变化情况，最终，基于变化并结合判定条件，识别是否发生灾害事件；而频率计算方法，首先，基于识别方法识别统计周期内的灾害发生次数，然后计算泥石流的发生频率。
1.基于遥感影像的泥石流灾害识别方法，其特征在于，包括以下步骤：
根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集；
根据目标泥石流流域所在区域的数字高程模型，提取目标泥石流流域的边界范围，并基于流域内坡度的变化将流域划分为子区域；
S2、根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数，然后，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数；
S3、基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集；
S4、基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件。
2.如权利要求1所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤S1中，根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集，包括：A1、从遥感影像数据集UL中，按自然年度分别提取待识别年度及其前后两年的遥感影像，并构建各年度的年度数据集ULY，下标Y表示年度；
A2、逐景判定年度数据集ULY中各影像数据Iij的所属月份，所述Iij表示第j年的第i景的遥感影像数据；然后，对连续三年的遥感影像数据，根据夏半年及冬半年的划分进行重构，构建待识别年度的夏半年的前后两个冬半年的数据集UDY，下标Y表示年度。
3.如权利要求2所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤A2中，在根据夏半年及冬半年的划分进行重构后，首先，将待录入冬半年数据集UDY的各影像数据Iij，逐个计算其云层覆盖率，然后，仅将满足云层覆盖率条件的影像数据录入数据集，完成两个冬半年数据集UDY的构建。
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、根据影像数据Iij的光谱信息，逐像素的进行云层像素判定，并将满足如下组合条件的像素认定为云层像素：其中， 表示第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、根据影像数据Iij的云层像素的数量，按如下公式计算其云层覆盖率：其中，Cij为第j年第i景影像数据的云层覆盖率，Nij‑Cloud为影像数据Iij中云层像素的数量，Nij‑All为影像数据Iij中的总像素数量。
5.如权利要求1所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤S1中，根据目标泥石流流域所在区域的数字高程模型，采用ArcGIS软件，提取目标泥石流流域的边界范围。
6.如权利要求1所述的基于遥感影像的泥石流灾害识别方法，其特征在于，
步骤S1中，基于流域内坡度的变化将流域划分为子区域，包括：
B1、根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集；
B2、遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点；
B3、基于分区突变点，计算子区域的分界线，并将边界范围内的目标泥石流流域划分为子区域；
B4、通过匹配数字高程模型和卫星影像数据的位置信息，提取卫星影像数据中与目标泥石流流域的边界范围、分界线及主河道对应像素的坐标。
7.如权利要求6所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤B1中，根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集，包括：B11、遍历边界范围内目标泥石流流域每个点位的高程值，并将最高高程点emax作为目标泥石流流域的主沟道起始点D1；
B12、将主沟道起始点D1的空间位置及高程值录入主沟道点集，并将主沟道起始点D1作为中心点DC的初始输入；
B13、根据输入的中心点DC，遍历其周围的八个点位，将其中高程值最低的点位识别为主沟道点Dk，将Dk的空间位置及高程值录入主沟道点集，下标k为主沟道点的序号；
B14、判定是否到达目标泥石流流域的边界，若是，则完成流域的主沟道的识别；否则，将主沟道点Dk作为新的中心点DC，返回步骤B13。
8.如权利要求6所述的基于遥感影像的泥石流灾害识别方法，其特征在于，目标泥石流流域的子区域包括形成区、流通区和堆积区，步骤B2中，遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点，包括：B21、从主沟起始点D1开始，依序计算主沟道点集中相邻两点间的纵坡比降：其中，Sk为主沟道点Dk与Dk+1两点间的纵坡比降，∆hk为主沟道点Dk与Dk+1两点间的高程差值，∆lk为主沟道点Dk与Dk+1两点间的水平距离，下标k为主沟道点的序号；
B22、将主沟道进行分段，并基于各分段所包含主沟道点的纵坡比降，计算分段的平均纵坡比降：其中， 为第q个分段的平均纵坡比降，Sqs为第q个分段所包含的第s个主沟道点的纵坡比降，r为分段所包含主沟道点的数量；
B23、从主沟起始点D1开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于或等于50%时，将连续多个分段中的首个分段作为形成区和流通区的划界分段；然后，从形成区和流通区的划界分段开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于1%时，将连续多个分段中的首个分段作为流通区和堆积区的划界分段；
B24、基于形成区和流通区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A1；基于流通区和堆积区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A2，以Dmin‑A1和Dmin‑A2作为分区突变点。
9.如权利要求8所述的基于遥感影像的泥石流灾害识别方法，其特征在于，
步骤B24中，从形成区和流通区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A1和纵坡比降最大的主沟道点Dmax‑A1；从流通区和堆积区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A2和纵坡比降最大的主沟道点Dmax‑A2；且提取的主沟道点Dmin‑A1、Dmax‑A1、Dmin‑A2和Dmax‑A2的排序应满足：其中，D1为主沟道的起始点，DK为主沟道的终点；
以Dmin‑A1和Dmin‑A2作为分区突变点，并将Dmax‑A1和Dmax‑A2作为辅助计算点；
在步骤B3中，在数字高程模型的水平投影内，将通过点Dmin‑A1且垂直于Dmin‑A1和Dmax‑A1连线的直线作为形成区和流通区的分界线；将通过点Dmin‑A2且垂直于Dmin‑A2和Dmax‑A2连线的直线作为流通区和堆积区的分界线。
10.如权利要求8或9任一项所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤B22中，将主沟道按长度80 120米进行分段，步骤B23中，所述连续多个分段包括：连续~三个分段。
11.如权利要求1 9任一项所述的基于遥感影像的泥石流灾害识别方法，其特征在于，~步骤S2，首先，按下式，根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数：其中， 和 为第j个冬半年的第i景影像数据的第m个像素在0.63 0.69μm波段和~
然后，按下式，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数：其中，j表示在前的冬半年， j+1表示在后的冬半年。
12.如权利要求11所述的基于遥感影像的泥石流灾害识别方法，其特征在于，步骤S3中，基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集，包括：S31、基于步骤S2中，由前后两个冬半年的第一景影像数据，计算获得的归一化植被指数NDVIm1j+1和物源扰动变化检测指数Indexm1j+1，按如下条件逐像素进行判定，并将满足条件的像素录入新增物源待定像素数据集PPD：其中，ThredNDVI为归一化植被指数的判定阈值，ThredIndex为物源扰动变化检测指数的判定阈值；
S32、基于步骤S2中，由前后两个冬半年的第一景以外的其他景影像数据，计算获得的归一化植被指数NDVImkj+1和物源扰动变化检测指数Indexmkj+1，分别将各景对应数据中与数据集PPD各像素相同坐标的像素的数据，按如下条件逐景逐像素进行判定：其中，k∈(1,M],M=min{Nj,Nj+1}，Nj为第j个冬半年影像数据的景数；
根据在各景中均满足条件的像素的坐标，对数据集PPD中的像素进行筛选，形成新增物源像素数据集；
S33、根据新增物源像素数据集所包含各像素的坐标，按目标泥石流流域各子区域间的分界线进性重构，获得各子区域的新增物源像素数据集。
13.如权利要求1 9任一项所述的基于遥感影像的泥石流灾害识别方法，其特征在于，~目标泥石流流域的子区域包括形成区、流通区和堆积区；
步骤S4中，基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件，包括：S41、将流通区新增物源像素数据集PLT中各像素与主河道的最小距离，分别与预设的流通区空间距离阈值ThredLT进行比较，若小于阈值ThredLT，则判定该像素为流通区泥石流活动有效堆积物源像素；然后，基于判定结果，构建流通区的有效堆积物源像素数据集PLT‑final；
将堆积区新增物源像素数据集PDJ中各像素与主河道的最小距离，分别与预设的堆积区空间距离阈值ThredDJ进行比较，若存在小于阈值ThredDJ的像素，则判定新增物源像素数据集PDJ所包含像素均为堆积区泥石流活动有效堆积物源像素，否则，判定新增物源像素数据集PDJ所包含像素均为堆积区泥石流活动无效堆积物源像素；然后，基于判定结果，构建堆积区的有效堆积物源像素数据集PDJ‑final；
S42、若流通区的有效堆积物源像素数据集PLT‑final的像素数量大于或等于阈值ThredLT‑pixel，或堆积区的有效堆积物源像素数据集PDJ‑final的像素数量大于或等于阈值ThredDJ‑pixel，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
14.如权利要求13所述的基于遥感影像的泥石流灾害识别方法，其特征在于，在步骤S42完成后，还包括：S43、若形成区新增物源像素数据集PXC的像素数量大于或等于阈值ThredXC‑pixel，且流通区的有效堆积物源像素数据集PLT‑final或堆积区的有效堆积物源像素数据集PDJ‑final的像素数量大于或等于阈值ThredLT‑DJ，则，进一步统计待识别年度连续15日最大累计降雨量PAcc15，若待识别年度的PAcc15超过待识别年度前3年同期的平均水平，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
15.基于遥感影像的泥石流灾害频率计算方法，其特征在于，
首先，根据权利要求1 14任一项所述的基于遥感影像的泥石流灾害识别方法，对指定~的统计周期内的年度，逐年识别目标泥石流流域是否发生泥石流；
然后，按如下公式计算目标泥石流流域的灾害发生频率FdebrisFlow：其中，Nevent为统计周期内目标泥石流流域发生泥石流灾害事件的次数，Nyear为统计周期的年份数。
基于遥感影像的泥石流灾害识别及频率计算方法
技术领域
[0001] 本发明涉及基于摄影测量的地质灾害调查领域，具体涉及基于遥感影像的泥石流灾害识别及频率计算方法。
背景技术
[0002] 泥石流是我国山区常见的地质灾害类型之一，对山区经济发展和人民安全构成严重威胁。由于山区泥石流沟数量众多、分布广泛，且每条泥石流沟发生泥石流灾害的频次均有所差异，精准掌握泥石流灾害频率信息，对于科学评估和防治泥石流灾害具有重要意义。
[0003] 当前，泥石流灾害的发生频率信息主要依靠两种方式获取：一、现场调查；二、文献调查。
[0004] 其中，现场调查，由于泥石流分布范围广阔，且部分区域存在通行障碍，人员现场调查费时费力，仅能够解决一部分位于公路沿线临近的泥石流沟的灾害频率信息获取问题，此外，相关信息只是被参与现场调查的相关人员掌握，难于共享，不利于区域联防联控。
[0005] 文献调查，则是从相关报道、新闻等文献资料中进行查询，能够快速获得一些泥石流沟的灾害发生频次信息，但这些信息能够覆盖的泥石流沟数量有限，仅限于少量影响重大的灾害事件，且频次信息受限于截止至资料发表时间之前，后续的灾害频次信息难以持续更新，对广域范围内泥石流灾害评估和防治工作的信息贡献较为乏力。
[0006] 因此，目前，尚未出现针对广域范围内泥石流灾害频率信息计算的深入研究。
发明内容
[0007] 本发明所要解决的技术问题是：提出一种基于遥感影像的泥石流灾害识别及频率计算方法，解决传统方式获取泥石流灾害频率信息存在的效率低、覆盖面不足、频次信息难以更新的问题。
[0008] 本发明解决上述技术问题采用的技术方案是：
[0009] 基于遥感影像的泥石流灾害识别方法，包括以下步骤：
[0010] S1、数据准备：
[0011] 根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集；
[0012] 根据目标泥石流流域所在区域的数字高程模型，提取目标泥石流流域的边界范围，并基于流域内坡度的变化将流域划分为子区域；
[0013] S2、根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数，然后，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数；
[0014] S3、基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集；
[0015] S4、基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件。
[0016] 具体的，步骤S1中，根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集，包括：
[0017] A1、从遥感影像数据集UL中，按自然年度分别提取待识别年度及其前后两年的遥感影像，并构建各年度的年度数据集ULY，下标Y表示年度；
[0018] A2、逐景判定年度数据集ULY中各影像数据Iij的所属月份，所述Iij表示第j年的第i景的遥感影像数据；然后，对连续三年的遥感影像数据，根据夏半年及冬半年的划分进行重构，构建待识别年度的夏半年的前后两个冬半年的数据集UDY，下标Y表示年度。
[0019] 进一步的，步骤A2中，在根据夏半年及冬半年的划分进行重构后，首先，将待录入冬半年数据集UDY的各影像数据Iij，逐个计算其云层覆盖率，然后，仅将满足云层覆盖率条件的影像数据录入数据集，完成两个冬半年数据集UDY的构建。
[0020] 具体的，提取的遥感影像数据，其光谱波段包括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七个波段；
~ ~ ~ ~
[0021] 在步骤A2中，仅将满足云层覆盖率<1%的影像数据录入数据集，并按如下步骤计算云层覆盖率：
[0022] A21、根据影像数据Iij的光谱信息，逐像素的进行云层像素判定，并将满足如下组合条件的像素认定为云层像素：
[0023]
[0024] 其中， 表示第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七个波段；
~ ~ ~
[0025] A22、根据影像数据Iij的云层像素的数量，按如下公式计算其云层覆盖率：
[0026]
[0027] 其中，Cij为第j年第i景影像数据的云层覆盖率，Nij‑Cloud为影像数据Iij中云层像素的数量，Nij‑All为影像数据Iij中的总像素数量。
[0028] 具体的，步骤S1中，根据目标泥石流流域所在区域的数字高程模型，采用ArcGIS软件，提取目标泥石流流域的边界范围。
[0029] 具体的，步骤S1中，基于流域内坡度的变化将流域划分为子区域，包括：
[0030] B1、根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集；
[0031] B2、遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点；
[0032] B3、基于分区突变点，计算子区域的分界线，并将边界范围内的目标泥石流流域划分为子区域；
[0033] B4、通过匹配数字高程模型和卫星影像数据的位置信息，提取卫星影像数据中与目标泥石流流域的边界范围、分界线及主河道对应像素的坐标。
[0034] 最优的，步骤B1中，根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集，包括：
[0035] B11、遍历边界范围内目标泥石流流域每个点位的高程值，并将最高高程点emax作为目标泥石流流域的主沟道起始点D1；
[0036] B12、将主沟道起始点D1的空间位置及高程值录入主沟道点集，并将主沟道起始点D1作为中心点DC的初始输入；
[0037] B13、根据输入的中心点DC，遍历其周围的八个点位，将其中高程值最低的点位识别为主沟道点Dk，将Dk的空间位置及高程值录入主沟道点集，下标k为主沟道点的序号；
[0038] B14、判定是否到达目标泥石流流域的边界，若是，则完成流域的主沟道的识别；否则，将主沟道点Dk作为新的中心点DC，返回步骤B13。
[0039] 最优的，目标泥石流流域的子区域包括形成区、流通区和堆积区，步骤B2中，遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点，包括：
[0040] B21、从主沟起始点D1开始，依序计算主沟道点集中相邻两点间的纵坡比降：
[0041]
[0042] 其中，Sk为主沟道点Dk与Dk+1两点间的纵坡比降，∆hk为主沟道点Dk与Dk+1两点间的高程差值，∆lk为主沟道点Dk与Dk+1两点间的水平距离，下标k为主沟道点的序号；
[0043] B22、将主沟道进行分段，并基于各分段所包含主沟道点的纵坡比降，计算分段的平均纵坡比降：
[0044]
[0045] 其中， 为第q个分段的平均纵坡比降，Sqs为第q个分段所包含的第s个主沟道点的纵坡比降，r为分段所包含主沟道点的数量；
[0046] B23、从主沟起始点D1开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于或等于50%时，将连续多个分段中的首个分段作为形成区和流通区的划界分段；然后，从形成区和流通区的划界分段开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于1%时，将连续多个分段中的首个分段作为流通区和堆积区的划界分段；
[0047] B24、基于形成区和流通区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A1；基于流通区和堆积区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A2，以Dmin‑A1和Dmin‑A2作为分区突变点。
[0048] 最优的，步骤B24中，从形成区和流通区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A1和纵坡比降最大的主沟道点Dmax‑A1；从流通区和堆积区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A2和纵坡比降最大的主沟道点Dmax‑A2；且提取的主沟道点Dmin‑A1、Dmax‑A1、Dmin‑A2和Dmax‑A2的排序应满足：
[0049]
[0050] 其中，D1为主沟道的起始点，DK为主沟道的终点；
[0051] 以Dmin‑A1和Dmin‑A2作为分区突变点，并将Dmax‑A1和Dmax‑A2作为辅助计算点；
[0052] 在步骤B3中，在数字高程模型的水平投影内，将通过点Dmin‑A1且垂直于Dmin‑A1和Dmax‑A1连线的直线作为形成区和流通区的分界线；将通过点Dmin‑A2且垂直于Dmin‑A2和Dmax‑A2连线的直线作为流通区和堆积区的分界线。
[0053] 最优的，步骤B22中，将主沟道按长度80 120米进行分段，步骤B23中，所述连续多~个分段包括：连续三个分段。
[0054] 具体的，步骤S2，首先，按下式，根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数：
[0055]
[0056] 其中， 和 为第j个冬半年的第i景影像数据的第m个像素在0.63 0.69μm波~段和0.76 0.90μm波段的TOA值；
~
[0057] 然后，按下式，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数：
[0058]
[0059] 其中，j表示在前的冬半年， j+1表示在后的冬半年。
[0060] 具体的，步骤S3中，基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集，包括：
[0061] S31、基于步骤S2中，由前后两个冬半年的第一景影像数据，计算获得的归一化植被指数NDVIm1j+1和物源扰动变化检测指数Indexm1j+1，按如下条件逐像素进行判定，并将满足条件的像素录入新增物源待定像素数据集PPD：
[0062]
[0063] 其中，ThredNDVI为归一化植被指数的判定阈值，ThredIndex为物源扰动变化检测指数的判定阈值；
[0064] S32、基于步骤S2中，由前后两个冬半年的第一景以外的其他景影像数据，计算获得的归一化植被指数NDVImkj+1和物源扰动变化检测指数Indexmkj+1，分别将各景对应数据中与数据集PPD各像素相同坐标的像素的数据，按如下条件逐景逐像素进行判定：
[0065]
[0066] 其中，k∈(1,M],M=min{Nj,Nj+1}，Nj为第j个冬半年影像数据的景数；
[0067] 根据在各景中均满足条件的像素的坐标，对数据集PPD中的像素进行筛选，形成新增物源像素数据集；
[0068] S33、根据新增物源像素数据集所包含各像素的坐标，按目标泥石流流域各子区域间的分界线进性重构，获得各子区域的新增物源像素数据集。
[0069] 具体的，目标泥石流流域的子区域包括形成区、流通区和堆积区；
[0070] 步骤S4中，基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件，包括：
[0071] S41、将流通区新增物源像素数据集PLT中各像素与主河道的最小距离，分别与预设的流通区空间距离阈值ThredLT进行比较，若小于阈值ThredLT，则判定该像素为流通区泥石流活动有效堆积物源像素；然后，基于判定结果，构建流通区的有效堆积物源像素数据集PLT‑final；
[0072] 将堆积区新增物源像素数据集PDJ中各像素与主河道的最小距离，分别与预设的堆积区空间距离阈值ThredDJ进行比较，若存在小于阈值ThredDJ的像素，则判定新增物源像素数据集PDJ所包含像素均为堆积区泥石流活动有效堆积物源像素，否则，判定新增物源像素数据集PDJ所包含像素均为堆积区泥石流活动无效堆积物源像素；然后，基于判定结果，构建堆积区的有效堆积物源像素数据集PDJ‑final；
[0073] S42、若流通区的有效堆积物源像素数据集PLT‑final的像素数量大于或等于阈值ThredLT‑pixel，或堆积区的有效堆积物源像素数据集PDJ‑final的像素数量大于或等于阈值ThredDJ‑pixel，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
[0074] 进一步的，在步骤S42完成后，还包括：
[0075] S43、若形成区新增物源像素数据集PXC的像素数量大于或等于阈值ThredXC‑pixel，且流通区的有效堆积物源像素数据集PLT‑final或堆积区的有效堆积物源像素数据集PDJ‑final的像素数量大于或等于阈值ThredLT‑DJ，则，进一步统计待识别年度连续15日最大累计降雨量PAcc15，若待识别年度的PAcc15超过待识别年度前3年同期的平均水平，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
[0076] 另外，本发明还提供了一种基于遥感影像的泥石流灾害频率计算方法，包括：
[0077] 首先，根据如上所述的基于遥感影像的泥石流灾害识别方法，对指定的统计周期内的年度，逐年识别目标泥石流流域是否发生泥石流；
[0078] 然后，按如下公式计算目标泥石流流域的灾害发生频率FdebrisFlow：
[0079]
[0080] 其中，Nevent为统计周期内目标泥石流流域发生泥石流灾害事件的次数，Nyear为统计周期的年份数。
[0081] 本发明的有益效果是：
[0082] 本发明基于遥感影像，利用像元光谱差异化特征计算，能够实现灾害发生前泥石流沟道物源区松散物源自动识别，及灾害发生后在泥石流沟道及沟口堆积物质的自动提取；结合地貌形态特征因素，通过时间序列对比计算，精准提取时间周期内堆积物变化情况，实现灾害发生、规模体量的量化判定，进而实现泥石流沟灾害发生频率的快速计算。该方法适用于广域范围的大量泥石流沟的灾害识别及频率信息提取，革新了传统的泥石流频率信息获取方式。
附图说明
[0083] 图1为本发明基于遥感影像的泥石流灾害识别方法的实现框架图。
具体实施方式
[0084] 本发明旨在提出一种适用于广域范围内的基于遥感影像的泥石流灾害识别方法，及基于其的泥石流灾害频率计算方法，以解决传统方式获取泥石流灾害频率信息存在的效率低、覆盖面不足、频次信息难以更新的问题。
[0085] 本发明的基于遥感影像的泥石流灾害识别方法，包括：
[0086] S1、数据准备：
[0087] 根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集；
[0088] 根据目标泥石流流域所在区域的数字高程模型，提取目标泥石流流域的边界范围，并基于流域内坡度的变化将流域划分为子区域；
[0089] S2、根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数，然后，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数；
[0090] S3、基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集；
[0091] S4、基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件。
[0092] 下面结合实施例，对本发明的方案作进一步的描述。
[0093] 实施例：
[0094] 如图1所示，本实施例的广域范围泥石流灾害的识别及频率计算方法，包括以下步骤：
[0095] S1、数据准备，其包括地形数据和遥感影像数据的准备。
[0096] A、遥感影像数据的准备：
[0097] 根据待识别年度，提取连续三年的目标泥石流流域的遥感影像，所述连续三年为待识别年度及其前后两年；基于连续三年的遥感影像，构建待识别年度的夏半年的前后两个冬半年的数据集。
[0098] 由于泥石流主要由降雨导致，通常发生在夏半年，因此，泥石流灾害的发生，可通过其前后两个冬半年的影像对比加以识别。在北半球，夏半年，通常指春分日到秋分日之间的这段时间；冬半年，通常指秋季10月经冬季到来年春季3月的这段时间，为了方便统一计算，在本实施例中，冬半年具体指当年11月、12月和来年的1月和2月这段时间。
[0099] 具体的，包括：
[0100] A1、从遥感影像数据集UL中，按自然年度分别提取待识别年度及其前后两年的遥感影像，并构建各年度的年度数据集ULY，下标Y表示年度；
[0101] A2、逐景判定年度数据集ULY中各影像数据Iij的所属月份，所述Iij表示第j年的第i景的遥感影像数据；然后，对连续三年的遥感影像数据，根据夏半年及冬半年的划分进行重构，构建待识别年度的夏半年的前后两个冬半年的数据集UDY，下标Y表示年度。
[0102] 如待识别年度为2003年，则提取2002、2003和2004连续三年的数据，并构建2002年11月 2003年2月以及2003年11月 2004年2月两个冬半年的数据集。上述的景是卫星遥感影~ ~
像数据的单位，1景即卫星拍摄的1张照片，同一区域卫星在某一天过境只拍摄1张照片，即1景影像，时间区间内，不同重访周期对应相应区域的影像数据景数=时间跨度总天数/重访周期天数。
[0103] 在高原及山区，冬半年通常属旱季，通常无云，但为了排除云层的干扰，进一步的，步骤A2中，在根据夏半年及冬半年的划分进行重构后，首先，将待录入冬半年数据集UDY的各影像数据Iij，逐个计算其云层覆盖率，然后，仅将满足云层覆盖率条件的影像数据录入数据集，完成两个冬半年数据集UDY的构建。云层覆盖率的计算可以采用现有任意云检测算法，其选择与遥感数据所包含的光谱信息有关。
[0104] 在本实施例中，提取的遥感影像数据，其光谱波段包括：波长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%的影像数据录入数据集。
[0105] 具体的，在步骤A2中，并按如下步骤计算云层覆盖率：
[0106] A21、根据影像数据Iij的光谱信息，逐像素的进行云层像素判定，并将满足如下组合条件的像素认定为云层像素：
[0107]
[0108] 其中， 表示第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七个波段；
~ ~ ~
[0109] A22、根据影像数据Iij的云层像素的数量，按如下公式计算其云层覆盖率：
[0110]
[0111] 其中，Cij为第j年第i景影像数据的云层覆盖率，Nij‑Cloud为影像数据Iij中云层像素的数量，Nij‑All为影像数据Iij中的总像素数量。
[0112] B、地形数据的准备：
[0113] 首先，根据目标泥石流流域所在区域的数字高程模型，提取目标泥石流流域的边界范围，并基于流域内坡度的变化将流域划分为子区域。目标泥石流流域的边界范围可以采用现有任意的算法，在本实施例中，采用ArcGIS软件，提取目标泥石流流域的边界范围。
[0114] 然后，基于流域内坡度的变化将流域划分为子区域，具体包括：
[0115] B1、根据目标泥石流流域边界范围内的高程变化，识别边界范围内目标泥石流流域的主沟道，构建主沟道点集。
[0116] 流域中，一定是由高往低流动，因此，可以采用逐点比较的方式来识别主沟道，但其同时受到分辨率及微地形的影响，在分辨率较高的情况下，也可以采用多点平均高程或坡度的方式来克服微地形的影响，
[0117] 在本实施例中，数字高程模型的分辨率为30米，因此，采用逐点比较的方式，具体的，步骤B1包括：
[0118] B11、遍历边界范围内目标泥石流流域每个点位的高程值，并将最高高程点emax作为目标泥石流流域的主沟道起始点D1；
[0119] B12、将主沟道起始点D1的空间位置及高程值录入主沟道点集，并将主沟道起始点D1作为中心点DC的初始输入；
[0120] B13、根据输入的中心点DC，遍历其周围的八个点位，将其中高程值最低的点位识别为主沟道点Dk，将Dk的空间位置及高程值录入主沟道点集，下标k为主沟道点的序号；
[0121] B14、判定是否到达目标泥石流流域的边界，若是，则完成流域的主沟道的识别；否则，将主沟道点Dk作为新的中心点DC，返回步骤B13。
[0122] B2、遍历主沟道点集，计算主沟道点集各点间的纵坡比降，根据纵坡比降的变化，识别并提取分区突变点。
[0123] 子区域的划分可以根据需要进行划分，而要实现灾害发生前泥石流沟道物源区松散物源自动识别，及灾害发生后在泥石流沟道及沟口堆积物质自动提取，因此，最好的，将目标泥石流流域的子区域包括形成区、流通区和堆积区，具体的，步骤B2包括：
[0124] B21、从主沟起始点D1开始，依序计算主沟道点集中相邻两点间的纵坡比降：
[0125]
[0126] 其中，Sk为主沟道点Dk与Dk+1两点间的纵坡比降，∆hk为主沟道点Dk与Dk+1两点间的高程差值，∆lk为主沟道点Dk与Dk+1两点间的水平距离，下标k为主沟道点的序号；
[0127] B22、将主沟道进行分段，并基于各分段所包含主沟道点的纵坡比降，计算分段的平均纵坡比降：
[0128]
[0129] 其中， 为第q个分段的平均纵坡比降，Sqs为第q个分段所包含的第s个主沟道点的纵坡比降，r为分段所包含主沟道点的数量；
[0130] B23、从主沟起始点D1开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于或等于50%时，将连续多个分段中的首个分段作为形成区和流通区的划界分段；然后，从形成区和流通区的划界分段开始，依次判定各分段的平均纵坡比降，当首次出现连续多个分段的平均纵坡比降均小于1%时，将连续多个分段中的首个分段作为流通区和堆积区的划界分段；
[0131] B24、基于形成区和流通区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A1；基于流通区和堆积区的划界分段，识别并提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A2，以Dmin‑A1和Dmin‑A2作为分区突变点。
[0132] 上述的判定条件：平均纵坡比降均小于或等于50%、平均纵坡比降均小于1%，是根据川西地区已有的泥石流灾害样本统计获得的形成区、流通区和堆积区之间的分界地形判定阈值条件。当然，针对其他的区域，可以根据该区域的具体情况，如已有的泥石流灾害样本，进行设置；同时，若对流域子区域的划分采用其他方式，阈值条件应与之相匹配。
[0133] 上述的，将主沟道进行分段，同时基于连续多个分段判定，是为了避免微观地形，比如巨石或者断崖等的影响，而分段的长度以及连续多个分段的数量则与数字高程模型的分辨率有关。最优的，步骤B22中，将主沟道按长度80 120米进行分段，步骤B23中，所述连续~多个分段包括：连续三个分段。在实施例中，数字高程模型的分辨率为30米，每段主沟道包括3个点位，由于识别的主沟道点连线呈折线状，因此，3个点位，能覆盖约100米的范围。
[0134] B3、基于分区突变点，计算子区域的分界线，并将边界范围内的目标泥石流流域划分为子区域。分界线也即经过分区突变点的线，并基于该分界线将流域划分为子区域，其可根据需要进行设定，比如，按主沟道两侧坡地的坡度设置斜线进行分界。
[0135] 在本实施例中，为了方便计算，具体的，包括：
[0136] 步骤B24中，从形成区和流通区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A1和纵坡比降最大的主沟道点Dmax‑A1；从流通区和堆积区的划界分段中，提取分段所包含主沟道点中纵坡比降最小的主沟道点Dmin‑A2和纵坡比降最大的主沟道点Dmax‑A2；且提取的主沟道点Dmin‑A1、Dmax‑A1、Dmin‑A2和Dmax‑A2的排序应满足：
[0137]
[0138] 其中，D1为主沟道的起始点，DK为主沟道的终点；上述的(D1，Dmin‑A1］表示Dmax‑A1位于主沟起始点D1和Dmin‑A1之间，其中，圆括号表示不包括该点，方括号表示包括该点，其他区间同理。
[0139] 以Dmin‑A1和Dmin‑A2作为分区突变点，并将Dmax‑A1和Dmax‑A2作为辅助计算点；
[0140] 在步骤B3中，在数字高程模型的水平投影内，将通过点Dmin‑A1且垂直于Dmin‑A1和Dmax‑A1连线的直线作为形成区和流通区的分界线；将通过点Dmin‑A2且垂直于Dmin‑A2和Dmax‑A2连线的直线作为流通区和堆积区的分界线。
[0141] B4、通过匹配数字高程模型和卫星影像数据的位置信息，提取卫星影像数据中与目标泥石流流域的边界范围、分界线及主河道对应像素的坐标。
[0142] S2、首先，根据前后两个冬半年的数据集中遥感影像对应波段的TOA值，获取各遥感影像的像素级的归一化植被指数。
[0143] 归一化植被指数，英文缩写为NDVI，能反映出植物冠层的背景影响，如土壤、潮湿地面、雪、枯叶、粗糙度等信息，且与植被覆盖有关。因此，基于前后两个冬半年的归一化植被指数，能获得各像素的变化情况，也即后述的像素级的物源扰动变化检测指数。
[0144] 归一化植被指数的计算，如下式所述：
[0145]
[0146] 其中， 和 为第j个冬半年的第i景影像数据的第m个像素在0.63 0.69μm波~段和0.76 0.90μm波段的TOA值；
~
[0147] 然后，按下式，通过对前后两个冬半年遥感影像逐景逐像素对比的方式，根据各遥感影像的归一化植被指数计算获得像素级的物源扰动变化检测指数：
[0148]
[0149] 其中，j表示在前的冬半年， j+1表示在后的冬半年。
[0150] S3、基于像素级的物源扰动变化检测指数和归一化植被指数，识别目标泥石流流域的新增物源像素，并根据各子区域的分界，构建目标泥石流流域各子区域的新增物源像素数据集。
[0151] 为了降低计算量，在本实施例中，具体的，步骤S3包括：
[0152] S31、基于步骤S2中，由前后两个冬半年的第一景影像数据，计算获得的归一化植被指数NDVIm1j+1和物源扰动变化检测指数Indexm1j+1，按如下条件逐像素进行判定，并将满足条件的像素录入新增物源待定像素数据集PPD：
[0153]
[0154] 其中，ThredNDVI为归一化植被指数的判定阈值，ThredIndex为物源扰动变化检测指数的判定阈值；
[0155] S32、基于步骤S2中，由前后两个冬半年的第一景以外的其他景影像数据，计算获得的归一化植被指数NDVImkj+1和物源扰动变化检测指数Indexmkj+1，分别将各景对应数据中与数据集PPD各像素相同坐标的像素的数据，按如下条件逐景逐像素进行判定：
[0156]
[0157] 其中，k∈(1,M],M=min{Nj,Nj+1}，Nj为第j个冬半年影像数据的景数；
[0158] 根据在各景中均满足条件的像素的坐标，对数据集PPD中的像素进行筛选，形成新增物源像素数据集；
[0159] S33、根据新增物源像素数据集所包含各像素的坐标，按目标泥石流流域各子区域间的分界线进性重构，获得各子区域的新增物源像素数据集。
[0160] 在本实施例中，上述的阈值ThredNDVI和ThredIndex，根据川西地区已有灾害样本统计分析得到，ThredNDVI设定为0.05，ThredIndex设定为0.45。当然，针对其他的区域，可以根据该区域的具体情况，如已有的泥石流灾害样本，进行设置。
[0161] S4、基于各子区域的新增物源像素数据集，根据预设的灾害阈值判定条件，识别目标泥石流流域是否在待识别年度发生泥石流灾害事件。
[0162] 灾害阈值判定条件，主要根据已有灾害样本统计分析得到，并同子区域的划分有关。而要实现灾害发生前泥石流沟道物源区松散物源自动识别，也即对形成区物源变化的识别；及灾害发生后在泥石流沟道及沟口堆积物质自动提取，也即，流通区和堆积区物源变化的识别，因此，最好的，将目标泥石流流域的子区域包括形成区、流通区和堆积区。
[0163] 因此，具体的，步骤S4包括：
[0164] S41、将流通区新增物源像素数据集PLT中各像素与主河道的最小距离，分别与预设的流通区空间距离阈值ThredLT进行比较，若小于阈值ThredLT，则判定该像素为流通区泥石流活动有效堆积物源像素；然后，基于判定结果，构建流通区的有效堆积物源像素数据集PLT‑final；
[0165] 将堆积区新增物源像素数据集PDJ中各像素与主河道的最小距离，分别与预设的堆积区空间距离阈值ThredDJ进行比较，若存在小于阈值ThredDJ的像素，则判定新增物源像素数据集PDJ所包含像素均为堆积区泥石流活动有效堆积物源像素，否则，判定新增物源像素数据集PDJ所包含像素均为堆积区泥石流活动无效堆积物源像素；然后，基于判定结果，构建堆积区的有效堆积物源像素数据集PDJ‑final；
[0166] S42、若流通区的有效堆积物源像素数据集PLT‑final的像素数量大于或等于阈值ThredLT‑pixel，或堆积区的有效堆积物源像素数据集PDJ‑final的像素数量大于或等于阈值ThredDJ‑pixel，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
[0167] S43、若形成区新增物源像素数据集PXC的像素数量大于或等于阈值ThredXC‑pixel，且流通区的有效堆积物源像素数据集PLT‑final或堆积区的有效堆积物源像素数据集PDJ‑final的像素数量大于或等于阈值ThredLT‑DJ，则，进一步统计待识别年度连续15日最大累计降雨量PAcc15，若待识别年度的PAcc15超过待识别年度前3年同期的平均水平，则判定目标泥石流流域在待识别年度发生了泥石流灾害。
[0168] 同样的，根据对川西地区已有灾害样本的统计分析，ThredLT‑pixel设定为5，ThredDJ‑pixel设定为10； ThredXC‑pixel设定为5，ThredLT‑DJ设定为1。ThredLT设定为60米，ThredDJ设定为30米。而上述的各像素与主河道的最小距离，也即该像素与主河道对应像素间距离的最小值。
[0169] 连续15日最大累计降雨量PAcc15，也即一年中，以15日为单位，计算获得的累计降雨量的最大值，具体计算见下式：
[0170]
[0171] 上式中，下标mj表示第j年的m日。
[0172] 进一步的，基于上述的识别方法，能够实现基于遥感影像的泥石流灾害频率计算方法，包括：
[0173] 首先，根据前述的基于遥感影像的泥石流灾害识别方法，对指定的统计周期内的年度，逐年识别目标泥石流流域是否发生泥石流；
[0174] 然后，按如下公式计算目标泥石流流域的灾害发生频率FdebrisFlow：
[0175]
[0176] 其中，Nevent为统计周期内目标泥石流流域发生泥石流灾害事件的次数，Nyear为统计周期的年份数。
[0177] 尽管这里参照本发明的实施例对本发明进行了描述，上述实施例仅为本发明较佳的实施方式，本发明的实施方式并不受上述实施例的限制，应该理解，本领域技术人员可以设计出很多其他的修改和实施方式，这些修改和实施方式将落在本申请公开的原则范围和精神之内。
