CN106443701A - 基于时序水范围遥感影像的洪涝灾前预警方法 - Google Patents

Abstract

本发明公开了一种基于时序水范围遥感影像的洪涝灾前预警方法，包括如下步骤：1)输入多时相高分辨率地表反射率遥感影像并进行数据预处理；2)计算水体指数，并基于水体指数影像提取水体范围二值影像；3)提取待检测时期的水体特征及正常期水体特征；4)分别将待检测时期水体特征与正常期水体特征进行叠加对比分析，检测异常信息：5)分析跟踪检测到的异常信息，若观测区域原先本不应该有水体覆盖的地方被水体连续覆盖，说明该观测区发生洪涝灾害的风险增加，应予以警示。该方法从多个角度分析水体特征的异常，根据异常分析对观测区发生洪涝灾害的风险进行预警，从而有效解决洪涝灾害灾前监测预警的困难问题。

Description

基于时序水范围遥感影像的洪涝灾前预警方法

技术领域

本发明涉及洪涝灾害监测领域，特别是指一种基于时序水范围遥感影像的洪涝灾前预警方法。

背景技术

卫星遥感观测技术由于观测范围广、周期性重访等独特优势，已经应用到包括水利、环境保护、防灾减灾等多个领域。应用遥感技术进行洪涝灾害的监测，已经成为当前防灾减灾的研究重点。

目前，遥感技术在洪涝灾害监测领域的应用主要包括光谱遥感监测和微波遥感监测两种监测方法。对于光谱遥感监测，光学遥感影像的光谱信息量大，易于从中有效地提取水体信息。目前应用比较广泛的光学数据主要有TM、ETM、MSS、SPOT、CBERS、AVHRR、MODIS、CCD、HRV、HRG、MERSI等。EOS/MODIS数据具有时间分辨率高，重访周期短，可以获得泛洪区域优质遥感影像，并为收集多时相历史数据提供条件，同时该数据覆盖范围广，适合大面积洪水监测。MSS、TM、CBERS、SPOT等光学遥感数据中的空间分辨率比较高，一般被广泛用于洪灾发生前土地利用信息的提取，以进行洪涝监测的背景数据分析，以及洪涝发生前后水体变化信息的监测。

洪涝遥感监测评估离不开对灾前、灾中与灾后的连续性监测。洪涝灾害发生前主要进行常态化的动态监测，收集多时相的易洪涝多发区域的水体覆盖历史数据，对洪峰、降水过程进行预测预警；在洪涝灾害发生过程中，由于光谱遥感监测技术受到观测条件的限制，常常处于非常被动的地位，而微波遥感监测技术中，微波传感器能够有效穿透云雾，获取洪涝灾害期间地面的水情信息，微波遥感数据主要有SAR、ASAR、SSM/I、MRI等，同时由于水体对雷达波束的镜面反射，使得水体能够较好地从雷达影像中提取出来，因此，在灾中主要利用微波遥感技术对洪峰经过区域进行监控；洪涝灾害发生过后，主要是通过对灾前、灾中和灾后多时相遥感影像的水体覆盖范围的变化检测并结合社会经济数据进行灾情评估。

利用光谱传感器对洪涝易发区域进行长期动态监测，并结合遥感技术建立灾前预测模型，对洪水或河道的变化进行预测，已经成为卫星遥感灾前预测预警的重要组成部分。虽然不少学者在这方面仍然在进行不懈的努力，然而目前建立灾前监测预测模型，对洪涝灾害进行较为准确的灾前监测仍然存在较大的困难。因此在洪涝灾前预警应用方面还存在一定的难度，目前还没有形成较好的方法加以解决，尚处于探索阶段。

发明内容

本发明的目的在于提供一种能够有效进行洪涝灾前监测预警的基于时序水范围遥感影像的洪涝灾前预警方法。

为实现上述目的，本发明所提供的基于时序水范围遥感影像的洪涝灾前预警方法，包括如下步骤：

1)输入多时相高分辨率地表反射率遥感影像，结合历史洪涝灾害记录，剔除发生洪涝灾害的遥感影像，通过配准和相对辐射校正对输入的遥感影像进行数据预处理；多时相指的是研究区6年或6年以上多年同期或者年内某时段内的多景影像，高分辨率指的是遥感影像的空间分辨率要优于30m；

2)根据步骤1)中经过预处理的遥感影像，计算水体指数，并根据水体指数提取水体范围二值影像；所述水体指数为归一化水体指数NDWI、改进的归一化差异水体指数MNDWI或新型水体指数NWI中的一种；

3)根据步骤2)中的水体范围二值影像，提取待检测时期水体特征和正常期水体特征；其中，待检测时期水体特征包括待检测时期的水体最大覆盖范围、水体覆盖频率、水体覆盖时长和水体覆盖面积变化趋势规律；正常期水体特征包括正常期的水体最大覆盖范围、水体覆盖频率、水体覆盖时长和水体覆盖面积变化趋势规律；说明：本专利中，各项正常期水体特征的时间段应与待检测时期的时间段是对应的；各正常期水体特征在求解过程所涉及影像均不包括洪涝灾害时期的影像；

4)分别将待检测时期的水体覆盖范围、水体覆盖频率、水体覆盖时长、水体面积变化趋势与正常期的水体最大覆盖范围、水体覆盖频率、水体覆盖时长、水体覆盖面积变化趋势规律进行叠加对比分析(一一对比)，按下述方法检测异常信息：

4.1)将待检测时期的水体覆盖范围与正常期的水体最大覆盖范围进行叠置变化检测，若某像元点在待检测时期为有水，而在正常期无水，则该像元点作为异常点检出；

4.2)将研究区待检测时期内的水体覆盖时长、水体覆盖频率分别与正常期的水体覆盖频率、水体覆盖时长进行叠置统计分析，若待检测时期的值背离正常期的值且超出设定阈值(不同水体特征有不同阈值)，则该像元点作为异常点检出；

4.3)将待检测时期内的水体覆盖面积变化趋势与正常期的水体覆盖面积变化趋势规律进行对比分析，若待检测时期内的水体覆盖面积变化规律与正常期的变化规律相背离，则研究区的水体覆盖面积变化趋势异常；

5)分析跟踪步骤4)中检测的异常信息，若观测区域原先本不应该有水体覆盖的地方被水体连续覆盖，说明该观测区发生洪涝灾害的风险增加，应予以警示。

优选地，步骤2)中，采用归一化水体指数NDWI提取水体范围二值影像，包括如下步骤为：

A、根据步骤1)中经过预处理的遥感影像，利用NDWI计算公式，计算每个像素对应的NDWI值，获得由每个像素NDWI值组成的指数影像，所述NDWI计算公式为：

式中，p(Green)为绿波段反射率值，p(NIR)为近红外反射率值；

B、对于生成的NDWI指数影像，选取适当的阈值α进行二值分割，阈值α的取值区间为(-1,1)，进一步优选为[-0.5,0.5]，其具体合适的值还需要反复试验并结合实际的影像经过目视解译得出，以尽量准确地区分出水体像元为准；一般情况下，α的经验值是0。指数影像中像元值大于等于阈值α的像元被提取为水体像元，其他像元为非水体像元，水体像元取值为1，非水体像元取值为0，得到的二值影像即为水体范围二值影像。

优选地，步骤3)中，所述正常期的水体最大覆盖范围NA(Normal Area)按如下方式确定：

A、按下式对研究区多年同期的水体范围二值影像进行叠加分析：

式中，A_k(i,j)为研究区多年同期的水体范围二值影像，n₁为年数，i，j为对应图像的行列号；

B、设置适当的阈值β，β的取值范围为0～20％，提取出正常期水体最大覆盖范围影像NA(i,j)；如果B(i,j)的值大于β，则NA(i,j)＝1，否则NA(i,j)＝0。待检测时期的水体覆盖范围可以直接从当前水体覆盖范围二值影像获得。

本发明中，类似NA(i,j)这种描述均表示图像NA中行列号为i,j的像素点处的取值，NT、NF等依此类推，不再特别说明。

优选地，步骤3)中，所述正常期的水体覆盖时长(Normal Time)按如下方式确定：

A、提取水体覆盖时长(Time)

按照如下公式提取水体覆盖时长影像T(i,j)：

式中，A_k(i,j)表示的是研究区在年内观测时段内的水体范围二值影像，n₂为该时段影像个数，i，j为像元点对应的行列号，m为每一景影像代表的时间尺度；

当前待检测时期的水体覆盖时长也根据此算法计算；

B、提取正常期水体覆盖时长(Normal Time)

按照如下公式提取正常期水体覆盖时长影像NT(i,j)：

式中，T_k(i,j)为研究区多年同期的水体覆盖时长影像，n₁为年数，i，j为像元点对应的行列号。

优选地，步骤3)中，所述正常期的水体覆盖频率按如下方式确定：

A、提取水体覆盖频率(Frequent)

按照如下公式提取水体覆盖频率影像F(i,j)：

式中，A_k(i,j)为年内观测时段内的有时序(按时间序列排列)水体范围二值影像，A₁(i,j)为像元点(i,j)在年内观测时段的开始值，i，j为像元点对应的行列号，n₂为该观测时段内的影像个数；将每个像元点在观测时段内的像元值变化用一维序列X(i,j)表示，对于像元点(i,j)，X(i,j)＝{该像元点在影像1像元值，该像元点在影像2像元值，…，该像元点在影像n₂像元值}，各影像按时间顺序排列；依次向后检索一维序列X(i,j)，每次检索2个连续元素，步长为1；若检索到的两个元素为{0,1}，则f(0,1)＝1，否则f(0,1)＝0；

并根据此算法，计算待检测时期的水体覆盖时长。

B、提取正常期水体覆盖频率(Normal Frequent)

按照如下公式提取正常期水体覆盖频率影像NF(i,j)：

式中，F_k(i,j)为多年同期水体覆盖频率数据影像，n₁为年数，i，j为像元点对应的行列号。

优选地，步骤4.1)中，所述的水体覆盖范围异常(Abnormal Area)，按照如下公式提取水体覆盖范围异常影像AbA(i,j)：

AbA(i,j)＝A(i,j)-(A(i,j)∪NA(i,j))

式中，A(i,j)，NA(i,j)分别为待检测时期的水体范围二值影像、正常期的水体最大覆盖范围影像，i，j为对应图像的行列号，若AbA(i,j)为1，则为异常点。

优选地，步骤4.2)中，所述水体覆盖时长异常(Abnormal Time)，按照如下公式提取水体覆盖时长异常影像AbT(i,j)：

式中，T(i,j)、NT(i,j)分别为待检测时期的水体覆盖时长影像、正常期水体覆盖时长影像。i，j为对应图像的行列号。若阈值的范围在20％～100％之间，则该像元点作为异常检出；取值越大，表示水体覆盖时长的偏离程度越大，异常程度越大。

优选地，步骤4.2)中，所述水体覆盖频率异常(Abnormal Frequent)，按照如下公式提取水体覆盖频率异常影像AbF(i,j)：

式中，F(i,j)、NF(i,j)分别为待检测时期的水体覆盖频率影像、正常期水体覆盖频率影像。i，j为对应图像的行列号。若AbF(i,j)>γ，阈值γ的范围在20％～100％之间，则该像元点作为异常检出；γ取值越大，表示水体覆盖频率的偏离程度越大，异常程度越大。

优选地，步骤3)中，利用动态谐波回归模型DHR，拟合水体面积的正常季节周期规律，得到正常期的水体覆盖面积变化趋势规律；

步骤4.3)中，所述的水体覆盖面积动态趋势异常，按如下公式判断是否为异常信息：

式中，表示t_i时间节点真值与拟合值的差值，表示待检测的t_i时间节点的水体覆盖面积值，表示采用动态谐波回归模型DHR模拟得到的t_i时间节点的正常期水体覆盖面积值，t_i时间节点表示在观测时段内的任一一个时间点，若则为异常。

优选地，该方法还包括如下步骤：采用动态谐波回归模型(DHR)，基于步骤2)中提取的水体范围二值影像进行时空分析，以获得研究区水体覆盖面积的时空变化规律，包括水体面积的长期变化趋势，水体面积变化的季节周期规律以及拟合出的水体面积多年的正常变化趋势规律。

本发明的有益效果是：本发明通过光谱传感器对洪涝易发区域进行长时间动态连续监测，并获取连续的遥感影像数据，建立相关模型，获取水体覆盖范围、水体覆盖频率、水体覆盖时长以及水体覆盖面积的变化规律这几个水体特征的正常期变化规律，从多个角度来分析水体特征的异常，根据异常分析对观测区发生洪涝灾害的风险进行预警，从而有效解决洪涝灾害灾前监测预警的困难问题。

附图说明

图1为实施例中所提供的洪涝灾前预警方法的总体流程图。

图2为实施例中多时相多角度水体异常检测技术流程图。

图3为实施例中水体覆盖范围异常检测的结果图。

图4为实施例中水体覆盖时长异常检测的结果图。

图5为实施例中水体覆盖频率异常检测的结果图。

图6为实施例中水体覆盖面积趋势异常检测及时空分析结果图。

具体实施方式

下面结合附图和具体实施例对本发明作进一步的详细说明。

如图1～2所示，本发明所提供的基于时序水范围遥感影像的洪涝灾前预警方法，包括如下步骤：

1)输入多时相高分辨率地表反射率遥感影像，结合历史灾前记录，剔除发生洪涝灾害的数据，通过配准和相对辐射校正对输入的遥感影像进行数据预处理。

2)根据步骤1)获得的地表反射率影像计算水体指数，采用归一化水体指数NDWI提取水体范围，包括如下步骤为：

A、提取水体范围的方法有很多，现有的遥感影像水体提取方法主要分为以下四类：特征分割方法、模式分类方法、轮廓线检测方法以及多源数据综合分析方法。特征指数法是通过多光谱波段运算获取反映水体与背景反差的指数，进而利用直方图阈值分割得到提取结果，典型的水体特征指数如归一化水体指数(NDWI)、改进的归一化差异水体指数(MNDWI)、新型水体指数(NWI)等，其模型简单，并且有不错的提取结果。在这里我们使用NDWI阈值分割法来提取水体范围，根据输入的多时相的地表反射率遥感影像，利用NDWI计算公式，计算每个像素对应的NDWI值。

B、NDWI的计算表达式为：

其中，p(Green)为绿波段反射率值，p(NIR)为近红外反射率值。

对于生成的NDWI指数影像，选取适当的阈值α进行二值分割。选取适当的阈值α进行二值分割，阈值α的取值区间为(-1,1)，进一步优选为[-0.5,0.5]，其具体合适的值还需要反复试验并结合实际的影像经过目视解译得出，以尽量准确地区分出水体像元为准；一般情况下，α的经验值是0。指数影像中像元值大于等于阈值α的像元被提取为水体像元，其他像元为非水体像元，水体像元取值为1，非水体像元取值为0，得到的二值影像即为水体覆盖范围提取结果。

实施例中，根据鄱阳湖流域2009-2014年多年4-12月HJ1A卫星数据影像，4-12月每月一景，基于NDWI计算公式，计算这些影像的NDWI值；并根据NDWI影像，提取了研究区的水体范围。

3)根据步骤2)中的水体范围数据，按下述方法提取待检测时期水体特征和正常期水体特征：

3.1)根据步骤2)中的水体范围数据，提取水体覆盖频率、水体覆盖时长，并获取待检测时期的水体特征，包括待检测时期水体最大覆盖范围、待检测时期水体覆盖频率、待检测时期水体覆盖时长和待检测时期水体覆盖面积变化趋势规律；

3.2)根据步骤2)中的水体范围和步骤3.1)中的水体覆盖频率和水体覆盖时长，提取正常期水体特征，包括正常期水体最大覆盖范围、正常期水体覆盖频率、正常期水体覆盖时长和正常期水体覆盖面积变化趋势规律；

其中，正常期水体最大覆盖范围(Normal Area)是指根据研究区多年同期的水体覆盖范围影像，通过叠加分析，提取的在该时期研究区被水体覆盖的最大正常范围。

按如下方式确定正常期水体最大覆盖范围影像NA(i,j)：

A、按下式对研究区多年同期(多年的含义是6年或6年以上，同期的含义是指年内某日的一景影像或某旬或月内的一景多天合成影像)的水体范围影像进行叠加分析：

式中，A_k(i,j)为研究区多年同期的水体覆盖范围二值影像，n₁为年数，i，j为对应图像的行列号，计算时按照影像的时间顺序依次相加。

B、设置适当的阈值β，β的取值范围为0～20％，提取出研究区被水体覆盖的最大正常范围NA(i,j)，如果B(i,j)的值大于β，则NA(i,j)＝1，否则，NA(i,j)＝0。

阈值β的选取是为了剔除掉某像元点由于云覆盖或偶然降水导致的水体像元。

实施例中，基于HJ1A卫星数据获得的鄱阳湖流域2009-2014年多年6月的水体覆盖范围数据影像，根据历史洪涝灾害记录，2010年6月鄱阳湖流域发生了洪涝灾害，所以剔除了2010年的数据。利用步骤3.2)中正常期水体最大覆盖范围的算法，将其他年份的数据进行叠加分析，阈值a取为10，提取鄱阳湖6月的正常期最大水体覆盖范围。如图3所示，左图表示的是鄱阳湖在6月份未发生洪涝灾害情况时的最大覆盖面积；中间表示的是鄱阳湖2015年6月的水体覆盖范围。

步骤3.1)中，水体覆盖时长(Time)是指研究区在年内某观测时段内被水体覆盖的总时长。

按如下公式得到水体覆盖时长影像T(i,j)：

式中，A_k(i,j)表示的是研究区在年内该观测时段内的时间序列水体覆盖范围二值影像，n₂为该时段影像个数，i，j为像元点对应的行列号，m为每一景影像代表的时间尺度，计算时按照影像的时间顺序依次相加。

步骤3.2)中，正常期的水体覆盖时长(Normal Time)是指根据研究区多年同期的水体覆盖时长影像，通过叠加分析，计算研究区各像元点在该观测时段内被水覆盖的平均时长。

按如下公式得到正常期水体覆盖时长影像NT(i,j)：

式中，T_k(i,j)为研究区多年同期水体覆盖时长影像，n₁为年数，i，j为像元点对应的行列号。

实施例中，基于HJ1A卫星数据所获得的鄱阳湖流域2009-2014年4-12月的水体覆盖范围影像，4-12观测时间段内一共9景影像，每月一景，每一景影像代表的时间尺度为30天。根据步骤3.1)中的计算水体覆盖时长的算法，分别计算2009-2014年每年4-12月这一观测时段的水体覆盖时长；再利用步骤3.2)中的正常期水体覆盖时长算法，根据2009-2014年6年的水体覆盖时长影像计算鄱阳湖4-12月的正常期水体覆盖时长。如图4所示，图4中的左图表示的是在4-12月这个时间段，鄱阳湖流域在未发生洪涝灾害时，被水体覆盖的正常时长，为水体覆盖时长异常检测提供正常期的参照基础；图4中的中间图表示的是2015年4-12月的水体覆盖时长，为待检测时期的水体覆盖时长。

步骤3.1)中，水体覆盖频率(Frequent)是指研究区在年内观测时段内被水体覆盖变动的次数或频次。

按如下公式得到年内该观测时段的水体覆盖频率影像F(i,j)：

式中，年内某观测时段内有时序水体范围二值影像A_k(i,j)，A₁(i,j)表示像元点(i,j)在年内观测时段的开始值，i，j为像元点对应的行列号，n₂为该观测时段内的影像个数；将每个像元点在观测时段内的像元值变化用一维序列X(i,j)表示，对于像元点(i,j)，X(i,j)＝{该像元点在影像1像元值，该像元点在影像2像元值，…，该像元点在影像n₂像元值}；依次向后检索一维序列X(i,j)，每次检索2个连续元素，步长为1；若检索到的两个元素为{0,1}，则f(0,1)＝1，否则f(0,1)＝0。

步骤3.2)中所述的正常期的水体覆盖频率NF(Normal Frequent)是指根据研究区多年同期的水体覆盖频率影像，通过叠加分析，计算研究区各像元点在该观测时段内被水覆盖的平均频率。

按如下公式得到正常期水体覆盖频率影像NF(i,j)：

式中，F_k(i,j)为研究区多年同期水体覆盖频率数据影像，n₁为年数，i，j为像元点对应的行列号。

实施例中，基于HJ1A卫星数据所获得的鄱阳湖流域2009-2014年4-12月的水体覆盖范围影像，4-12观测时间段内一共9景影像，每月一景。根据步骤3.1)中计算水体覆盖频率的算法，分别计算2009-2014年每年年内的水体覆盖频率；再利用步骤3.2)中正常期水体覆盖频率的算法，根据2009-2014多年的水体覆盖频率影像计算正常期水体覆盖频率。如图5所示，图5中的左图表示的是鄱阳湖流域在4-12月这个时间段，未发生洪涝灾害时被水体覆盖的正常频率，为水体覆盖频率异常检测提供正常期的参照基础；图5中的中间图表示的是2015年4-12月这一观测时段内的水体覆盖频率，为待检测水体覆盖频率。

步骤3.2)中所述的正常期水体覆盖面积变化趋势规律是指研究区的水体覆盖面积在观测时段内的正常变化趋势规律。

按如下方式得到正常期水体覆盖面积变化趋势规律：

根据研究区长时间序列的水体覆盖范围二值影像，基于动态谐波回归模型(DHR)，拟合水体面积的正常变化趋势规律。

动态谐波回归模型(DHR)基于UC模型发展而来，UC模型有多种表现形式，DHR(Dynamic Harmonic Regession)模型可以认为是UC模型的一种特殊表现形式(Young等，1999)，它建立在UC模型的基础上。

对于一个时间序列，DHR模型主要用于拟合时间序列的3个成分，即T_t、S_t及e_t，表达为：

y_t＝T_t+S_t+e_t e_t～N(0,σ²)

式中，y_t表示原始时间序列在t时刻的值；T_t表示为原始序列在该时刻的趋势成分；S_t表示原始序列在该时刻的周期成分，其中S_t为季节周期成分；e_t为残差，通常为了简便，e_t被视作为高斯白噪声分布。本技术中的研究对象为水体面积，水体面积具有季节和周期变化规律，DHR模型的主要特点是对季节和周期成分进行拟合，因此它适用于有季节或周期变化成分的时间序列分析。

本技术中的DHR模型是以MATLAB平台的Captain Toolbox软件包为基础进行开发建模的，该软件包由Lancaster University的Young等设计并发展(http://www.es.lancs.ac.uk/cres/captain/)。根据开发的DHR模型功能，直接输入长时间序列水体覆盖范围影像，就可以拟合出研究区水体面积的正常变化趋势规律。

本实施例中，基于HJ1A卫星数据所获得的鄱阳湖流域2009-2014年4-12月的水体覆盖范围影像，4-12观测时间段内一共9景影像，每月一景。根据步骤3.2)中计算正常期水体覆盖面积变化趋势的算法，利用动态谐波回归模型(Dynamic Harmonic Regession)对研究区的水体覆盖面积进行分析。如图6所示有4条曲线，在本步骤中，需要的结果是拟合水体面积，表示的是研究区的真实水体面积在去除噪声的扰动后，拟合的水体覆盖面积的正常变化趋势规律，为水体面积趋势异常检测提供正常依据。

4)分别将待检测时期的水体覆盖范围、水体覆盖频率、水体覆盖时长和水体覆盖面积变化趋势与正常期水体最大覆盖范围、正常期水体覆盖频率、正常期水体覆盖时长以及正常期水体覆盖面积变化趋势规律进行叠加对比分析，检测异常信息。

4.1)水体覆盖范围异常检测是以传统异常检测方法理论为基础的异常检测方法，即在传统异常检测理论中，认为如数据与给定统计分布或模型的显著性差异超过某一特定数值或范围即为异常。在进行水体覆盖范围异常检测时，将当前待检测的水体范围影像与该时间的正常期最大水体覆盖范围进行叠置变化检测分析，提取出与正常期的水体分布不一致部分，即为异常的水体覆盖范围。

4.2)水体覆盖时长或频率异常检测是以基于传统异常检测方法理论为基础的异常检测方法，即在传统异常检测理论中，认为如数据与给定统计分布或模型的显著性差异超过某一特定数值或范围即为异常。在进行水体覆盖时长或频率异常检测时，将当前待检测的水体覆盖时长或频率影像与同期的正常期特征进行叠置阈值分析，若待检测时期的值超过正常期特征值一定阈值之外则认为是异常的。

4.3)水体覆盖面积动态趋势异常检测是以偏差异常检测理论中基于轮廓的方法为基础的一种异常检测方法。基于轮廓的异常检测方法，首先给定一个时间序列的正常轮廓，将新的时间点或时间段与轮廓进行比对来识别其异常与否。这里，轮廓实际上是对时间序列的历史数据进行平滑而得到的一个标准序列，即这类的方法的核心在于如何获得时间序列的正常轮廓。

步骤4.1)中所述的水体覆盖范围异常(Abnormal Area)，按照如下公式提取异常水体范围影像AbA(i,j)：

AbA(i,j)＝A(i,j)-(A(i,j)∪NA(i,j))

式中，A(i,j)，NA(i,j)分别为待检测时期的水体范围二值影像、正常期最大水体范围二值影像，i，j为对应图像的行列号，若AbA(i,j)为1，则为异常点。

实施例中，根据步骤3.2)获取的6月鄱阳湖流域正常期最大水体覆盖范围并结合2015年6月的水体范围影像，利用步骤4.1)中检测水体范围异常的算法，提取的异常水体范围如图3中右图所示。通过对研究区水体范围的异常检测，我们就能够知道待检测时期的水体范围是否超出了正常期的最大范围，若超出了，就应当及时预警。

步骤4.2)中所述的水体覆盖时长异常(Abnormal Time)，按照如下公式提取水体覆盖时长异常影像AbT(i,j)：

式中，T(i,j)、NT(i,j)分别为待检测水体覆盖时长影像、正常期水体覆盖时长影像。i，j为对应图像的行列号。若阈值的范围在20％-100％之间，则该像元点作为异常检出；取值越大，表示水体覆盖时长的偏离程度越大，异常程度越大。

实施例中，根据步骤3.2)获取的鄱阳湖流域4-12月的正常期水体覆盖时长并结合2015年4-12月的水体覆盖时长影像，利用步骤4.2)中检测水体覆盖时长异常的算法，阈值选择为20％，检测水体覆盖时长异常，并输出异常点的空间分布，结果如图4中右图所示。水体覆盖时长异常表示的是异常区域在研究时段内被水体淹没的时长超出了正常期的值，表明异常区域原先是没有水或者是较少被水覆盖的，应当及时预警。

步骤4.2)所述的水体覆盖频率异常(Abnormal Frequent)，按照如下公式提取水体覆盖频率异常影像AbF(i,j)：

式中，F(i,j)、NF(i,j)分别为待检测水体覆盖频率影像、正常期水体覆盖频率影像。i，j为对应图像的行列号。若AbF(i,j)>γ，阈值γ的范围在20％-100％之间，则该像元点作为异常检出；γ取值越大，表示水体覆盖频率的偏离程度越大，异常程度越大。

实施例中，根据步骤3.2)获取的鄱阳湖流域4-12月的正常期水体覆盖频率并结合2015年4-12月的水体覆盖频率影像，利用步骤4.2中检测水体覆盖频率异常的算法，阈值γ选择为20％，检测水体覆盖频率异常，并输出异常点的空间分布，结果分别如图5中右图所示。水体覆盖频率异常表示的是异常区域在研究时段内被水体频繁淹没，而异常表明这些异常区域的被水体淹没的风险提高，应当及时预警。

步骤4.3)中所述的水体覆盖面积动态趋势异常检测，按照如下方式提取异常信息：

式中，表示t_i时间节点真值与拟合值的差值，表示t_i时间节点的真实的水体面积值，表示t_i时间节点的拟合的正常期水体面积值，如果则为异常，否则为非异常。

实施例中，根据步骤3.2)的算法，所获取的鄱阳湖水体面积正常变化趋势规律，利用步骤4.3)中检测水体覆盖面积动态趋势异常的算法，将水体面积的真实值与拟合的正常趋势规律叠加分析，如图6所示，圆圈部分为真值高于拟合的正常规律的时间段，若异常持续的时间过长，则说明该区域被水体淹没的面积持续增大，发生洪涝的风险增加，应该及时预警。

步骤5，基于DHR模型，对研究区多年长时间序列的水体覆盖范围影像进行时空分析，探究研究区水体覆盖面积的时空变化规律。能够分析出研究区水体面积的多年长期变化趋势，水体面积年内的季节周期变化规律。

实施例中，如图6，趋势，即研究区在这段时间内水体面积总体变化的趋势；周期，即研究区的水体覆盖面积在年内的季节周期变化规律。从水体面积变化的季节周期规律可以看出，湖水面积会随着季节的变化呈现出丰水期和枯水期的交替变化，这与鄱阳湖是一个季节性吞吐湖是相符合的；另外从水体面积变化趋势线中可以看出，研究区在这段时间水体面积的整体变化趋势呈现轻微下降趋势，和近年来鄱阳湖水体面积逐渐萎缩的实际情况相符合。

Claims (10)

1.一种基于时序水范围遥感影像的洪涝灾前预警方法，其特征在于：包括如下步骤：

1)输入多时相高分辨率地表反射率遥感影像，结合历史洪涝灾害记录，剔除发生洪涝灾害的遥感影像，通过配准和相对辐射校正对输入的遥感影像进行数据预处理；

2)根据步骤1)中经过预处理的遥感影像，计算水体指数，并根据水体指数提取水体范围二值影像；所述水体指数为归一化水体指数NDWI、改进的归一化差异水体指数MNDWI或新型水体指数NWI；

3)根据步骤2)中的水体范围二值影像，提取待检测时期水体特征及正常期水体特征；其中：待检测时期水体特征包括待检测时期的水体覆盖范围、水体覆盖频率、水体覆盖时长和水体覆盖面积变化趋势规律；正常期水体特征包括正常期的水体最大覆盖范围、水体覆盖频率、水体覆盖时长和水体覆盖面积变化趋势规律；

4)分别将待检测时期的水体覆盖范围、水体覆盖频率、水体覆盖时长、水体面积变化趋势与正常期的水体最大覆盖范围、水体覆盖频率、水体覆盖时长、水体覆盖面积变化趋势规律进行叠加对比分析，按下述方法检测异常信息：

4.1)将待检测时期的水体覆盖范围与正常期的水体最大覆盖范围进行叠置变化检测，若某像元点在待检测时期为有水，而在正常期无水，则该像元点作为异常点检出；

4.2)将待检测时期的水体覆盖时长、水体覆盖频率分别与正常期的水体覆盖频率、水体覆盖时长进行叠置统计分析，若待检测时期的值背离正常期的值且超出设定阈值，则该像元点作为异常点检出；

4.3)将待检测时期的水体覆盖面积变化趋势与正常期的水体覆盖面积变化趋势规律进行对比分析，若待检测时期的水体覆盖面积变化规律与正常期的变化规律相背离，则水体覆盖面积变化趋势异常；

5)分析跟踪步骤4)中检测的异常信息，若观测区域原先本不应该有水体覆盖的地方被水体连续覆盖，说明该观测区发生洪涝灾害的风险增加，应予以警示。

2.根据权利要求1所述的基于时序水范围遥感影像的洪涝灾前预警方法，其特征在于：

步骤2)中，根据归一化水体指数NDWI提取水体范围二值影像，包括如下步骤：

A、根据步骤1)中经过预处理的遥感影像，利用NDWI计算公式，计算每个像素对应的NDWI值，获得由每个像素NDWI值组成的指数影像，所述NDWI计算公式为：

$NDWI=\frac{p\left(Green\right)-p\left(NIR\right)}{p\left(Green\right)+p\left(NIR\right)}$

式中，p(Green)为绿波段反射率值，p(NIR)为近红外反射率值；

B、对于生成的NDWI指数影像，选取阈值α进行二值分割，阈值α的取值区间为(-1,1)，指数影像中像元值大于等于阈值α的像元被提取为水体像元，其他像元为非水体像元，水体像元取值为1，非水体像元取值为0，得到的二值影像即为水体范围二值影像。

3.根据权利要求2所述的基于时序水范围遥感影像的洪涝灾前预警方法，其特征在于：

步骤3)中，所述正常期的水体最大覆盖范围按如下方式确定：

A、按下式对多年同期的水体范围二值影像进行叠加分析：

$B(i,j)=\frac{{\mathrm{\Σ}}_{k=1}^{n_1}{A}_k(i,j)}{n_1}$

式中，A_k(i,j)为多年同期的水体范围二值影像，n₁为年数，i，j为对应图像的行列号。

B、设置适当的阈值β，β的取值范围为0～20％，提取出正常期水体最大覆盖范围影像NA(i,j)；如果B(i,j)的值大于β，则NA(i,j)＝1，否则NA(i,j)＝0。

4.根据权利要求2所述的基于时序水范围遥感影像的洪涝灾前预警方法，其特征在于：

步骤3)中，所述正常期的水体覆盖时长按如下方式确定：

A、按照如下公式提取水体覆盖时长影像T(i,j)：

$T(i,j)=\left(\underset{k=1}{\overset{n_2}{\mathrm{\Σ}}}{A}_k\right(i,j\left)\right)*m$

式中，A_k(i,j)表示的是年内观测时段内的水体范围二值影像，n₂为该时段内影像个数，i，j为像元点对应的行列号，m为每一景影像代表的时间尺度；

B、按照如下公式提取正常期水体覆盖时长影像NT(i,j)：

$NT(i,j)=\frac{{\mathrm{\Σ}}_{k=1}^{n_1}{T}_k\left(i,j\right))}{n_1}$

式中，T_k(i,j)表示的是多年同期的水体覆盖时长影像，n₁为年数，i，j为像元点对应的行列号。

5.根据权利要求2所述的基于时序水范围遥感影像的洪涝灾前预警方法，其特征在于：

步骤3)中，所述正常期的水体覆盖频率按如下方式确定：

A、按照如下公式提取水体覆盖频率影像F(i,j)：

式中，A_k(i,j)为年内观测时段内的有时序水体范围二值影像，A₁(i,j)为像元点(i,j)在年内观测时段的开始值，i，j为像元点对应的行列号，n₂为该观测时段内的影像个数；将每个像元点在观测时段内的像元值变化用一维序列X(i,j)表示，对于像元点(i,j)，X(i,j)＝{该像元点在影像1像元值，该像元点在影像2像元值，…，该像元点在影像n₂像元值}；依次向后检索一维序列X(i,j)，每次检索2个连续元素，步长为1；若检索到的两个元素为{0,1}，则f(0,1)＝1，否则f(0,1)＝0；

B、按照如下公式提取正常期水体覆盖频率影像NF(i,j)：

$NF(i,j)=\frac{{\mathrm{\Σ}}_{k=1}^{n_1}\left(F_k\right(i,j\left)\right)}{n_1}$

式中，F_k(i,j)为多年同期水体覆盖频率影像，n₁为年数，i，j为像元点对应的行列号。

6.根据权利要求3所述的基于时序水范围遥感影像的洪涝灾前预警方法，其特征在于：

步骤4.1)中，按照如下公式提取水体覆盖范围异常影像AbA(i,j)：

AbA(i,j)＝A(i,j)-(A(i,j)∪NA(i,j))

式中，A(i,j)，NA(i,j)分别为待检测时期的水体覆盖范围影像、正常期水体最大覆盖范围影像，i，j为对应图像的行列号；若AbA(i,j)为1，则为异常点。

7.根据权利要求4所述的基于时序水范围遥感影像的洪涝灾前预警方法，其特征在于：

步骤4.2)中，按照如下公式提取水体覆盖时长异常影像AbT(i,j)：

$AbT(i,j)=\frac{|T(i,j)-NT(i,j)|}{\max \left\{NT(i,j)\right\}}$

式中，T(i,j)、NT(i,j)分别为待检测时期的水体覆盖时长影像、正常期水体覆盖时长影像，i，j为对应图像的行列号，若则该像元点作为异常检出，阈值的范围在20％～100％之间。

8.根据权利要求5所述的基于时序水范围遥感影像的洪涝灾前预警方法，其特征在于：

步骤4.2)中，按照如下公式提取水体覆盖频率异常影像AbF(i,j)：

$AbF(i,j)=\frac{|F(i,j)-NF(i,j)|}{\max \left\{NF(i,j)\right\}}$

式中，AbF(i,j)、AbF(i,j)分别为待检测时期的水体覆盖频率影像、正常期水体覆盖频率影像，i，j为对应图像的行列号，若AbF(i,j)>γ，则该像元点作为异常检出，阈值γ的范围在20％～100％之间。

9.根据权利要求2所述的基于时序水范围遥感影像的洪涝灾前预警方法，其特征在于：

步骤3)中，利用动态谐波回归模型，拟合水体面积的正常季节周期规律，得到正常期的水体覆盖面积变化趋势规律；

步骤4.3)中，进行水体覆盖面积动态趋势异常检测时，按如下公式判断是否为异常信息：

${\mathrm{\ΔS}}_{t_i}=S_{t_i}-{S_{t_i}}^{\′}$

式中，表示t_i时间节点真值与拟合值的差值，表示待检测的t_i时间节点的水体覆盖面积值，表示采用动态谐波回归模型模拟得到的t_i时间节点的正常期水体覆盖面积值，若则为异常。

10.根据权利要求1～9中任一项所述的基于时序水范围遥感影像的洪涝灾前预警方法，其特征在于：

采用动态谐波回归模型DHR，对步骤2)中提取的水体范围二值影像进行时空分析，以获得水体覆盖面积的时空变化规律，包括水体面积长期变化的趋势，水体面积变化的季节周期规律以及拟合出的水体面积多年的正常变化趋势规律。