一种基于SAR影像的洪涝风险区自动提取方法与流程

一种基于sar影像的洪涝风险区自动提取方法
技术领域
1.本发明属于洪涝监测技术领域，具体地说，发明了一种基于sar影像的洪涝风险区自动提取方法。该方法利用c波段后向散射系数(σ0)计算洪涝事件前后的两时相sar影像的水体概率，再利用两时相影像的水体概率增量表达降雨事件后新增洪涝淹没区域(洪涝风险区域)。该方法是一种简单、快捷、有效的洪涝监测方法，可以快速监测水体面积变化，自动提取洪涝风险区域。


背景技术：

2.全球气候变暖，异常气候现象频繁发生，洪涝灾害渐趋频繁，对人类的生命健康、生存环境造成的巨大影响。尤其我国南方地区，多处于热带和亚热带季风气候，强降雨和上游的冰雪消融现象对该区域的造成了更为严重的洪涝灾害。因此，有效快速的汛期水体变化检测，自动提取洪涝风险区域，十分重要。
3.在强降雨、阴雨天气下，光学影像成像条件很差，存在大量的数据缺失，制约了光学影像在汛期水体变化检测中的应用。合成孔径雷达(sar)不受阴雨影响，具有全天时全天候的对地观测能力，在水体范围提取和洪涝风险区域提取具有巨大潜力。目前，有大量在轨运行的sar卫星，如：sentinel
‑
1卫星、gf
‑
3卫星等，都是优越的sar影像数据源，为我国洪涝监测提供了重要数据支撑。
4.通常，基于遥感数据的洪涝风险区域提取是先提取两时相影像的水体部分，再进行叠置分析，分析洪灾前后水域变化。水体环境十分复杂的情况下，水体提取易出现误提和漏检等错分现象，在一定程度上，像元所属地物类型的不确定性是原因之一。另外，也有研究先对不同时期的影像后向散射进行运算，再利用阈值分割等方法对运算结果进行水体变化检测，进行灾后评估和洪涝风险区提取。由于水体环境的复杂情况，阈值分割也会具有一定的不确定性，而且，阈值设置也缺乏自动化。


技术实现要素：

5.本发明的目的是提供一种洪涝风险区域提取方法，直接基于洪涝事件前后的两时相sar影像进行洪涝监测，通过计算洪涝事件后的水体概率增量，构建基于水体概率的洪涝风险区域提取方法，可以快速进行洪涝监测。
6.为实现上述目的，本发明所采用的技术方案是：
7.一种基于sar影像的洪涝风险区域自动提取方法，其特征是在于，包括如下步骤：
8.步骤(1)，获取区域内洪涝事件发生前后的两时相影像，对两时相影像进行预处理，并将所有像元组成后向散射系数(σ0)集合σ1；
9.步骤(2)，对集合σ1绘制频率直方图，假设分布有水体的两时相sar影像的σ0集合符合混合高斯分布，并对σ0集合中的离群值过滤；
10.步骤(3)，对步骤(2)过滤后对应的两时相sar影像进行降低空间分辨率的重采样，生成σ0集合σ2，并假设降采样后的σ2呈混合高斯分布，且所形成的分布可用于降采样前集
合σ1的水体概率计算，对降采样后的σ0集合σ2进行k均值聚类，计算水体和非水体先验概率；
11.步骤(4)，利用最大期望(em)算法和步骤(3)的水体先验概率对将降采样后的集合σ2进行迭代计算分布参数，以提高迭代速度，然后对降采样前的每个像元的σ0值进行水体概率计算；
12.步骤(5)，以水体概率增量表达事件前后每个像元属于水体的概率的变化值，确定洪涝风险区与非风险区的水体概率增量边界值，构建基于两时相水体概率的洪涝风险区域提取方法。
13.在采用上述技术方案的基础上，本发明还可采用以下进一步的技术方案或对这些进一步的技术方案组合使用：
14.所述步骤(1)中具体包括以下子步骤：
15.步骤(11)，获取洪涝事件前后覆盖研究区的sar影像；
16.步骤(12)，利用sarscape软件对两时相sar影像进行配准、多时相滤波、地理编码、辐射定标、空间裁剪，获取σ0(单位：db)；
17.步骤(13)，将两时相影像所有像元进行组合，组成σ0集合σ1。
18.所述步骤(2)中对集合σ1绘制频率直方图，假设分布有水体的两时相sar影像的σ0集合符合混合高斯分布，并对离群值过滤，具体包括以下子步骤：
19.步骤(21)，对σ0集合绘制频率直方图，并假设其概率密度曲线符合混合高斯分布；
20.步骤(22)，即使已经利用地形数据和两时相特征对sar影像进行一系列预处理，但依然会存在较大的异常值影响后续参数估计和水体概率计算，因此，按数值从小到大排列，去除离群值。
21.所述步骤(3)中具体包括以下子步骤：
22.步骤(31)，为提高步骤(4)中的迭代计算速度，对对步骤(2)过滤后对应的两时相sar影像进行降采样生成σ0集合σ2，并假设集合σ2的概率密度函数也符合混合高斯分布，且基于集合σ2的概率密度函数可对步骤(1)的σ0计算水体概率。
23.步骤(32)，对σ0集合σ2进行k均值聚类，集合分为水体类的σ0集合和非水体类的σ0集合；
24.步骤(33)，计算水体和非水体先验概率；水体像元集合的像元个数占所有像元的个数比值为水体先验概率p0(w)，且p0(w)为p(w)的初值，非水体像元集合的像元个数占所有像元的个数比值为非水体先验概率，公式如下：
[0025][0026]
其中，w为水体像元集合，为非水体像元集合，p(w)与分别表示水体先验概率和非水体先验概率，且两者加和等于1。
[0027]
所述步骤(4)中具体包括以下子步骤：
[0028]
步骤(41)，利用混合高斯分布构建σ0边际分布的概率密度函数，密度函数如下：
[0029][0030]
其中，w为水体像元集合，为非水体像元集合，p(w)与分别表示水体先验概率和非水体先验概率，p(σ0)表示该后向散射系数对应的像元属于水体的概率；
[0031][0032][0033]
其中，μ
w
、s
w
、分别为水体像元σ0的均值、标准差、非水体像元σ0的均值、标准差；
[0034]
步骤(42)，利用基于最大似然估计的最大期望(em)算法对分布曲线的均值、标准差参数迭代计算，对p(w)进行迭代优化，公式如下：
[0035][0036]
其中，em(
·
)代表最大期望(em)算法，p(w)为由k均值聚类求得的水体先验概率，σ0为后向散射系数数据，μ
w
、s
w
、分别为水体像元σ0的均值、标准差、非水体像元σ0的均值、标准差。
[0037]
步骤(43)，将步骤(42)中优化的p(w)、p(σ0|w)和代入水体概率计算公式，构建水体概率计算公式，公式如下：
[0038]
p(w|σ0)＝p(w)p(σ0|w)/p(σ0)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0039]
步骤(44)，将降采样前的两时相影像的σ0代入公式，分别计算两时相影像的水体概率。
[0040]
所述步骤(5)中具体包括以下子步骤：
[0041]
步骤(51)，求取洪涝事件后前水体概率之差，其中，洪涝事件前的水体概率记为p1，洪涝事件后水体概率记为p2，用洪涝事件后的水体概率值减去洪涝事件前的水体概率值，记为p
21
；
[0042]
步骤(52)，风等会使后向散射较为粗糙，使得后向散射系数增大，但是通常情况下降雨事件前的水体在降雨事件后依然保持为水体，理想情况下p
21
不存在负值，因此取0为p
21
的关键阈值，利用水体概率进行洪涝风险区提取，并以概率增量的形式表达，提取规则如下：
[0043][0044]
其中，y大于0相对应的区域即为洪涝风险区域，y越接近1代表水体淹没状况变化越明显，y等于0的区域为非洪涝风险区域。
[0045]
步骤(53)，对y进一步进行阈值分割，设定0.2、0.4、0.6和0.8为阈值，0
‑
0.2对应风险一级，0.2
‑
0.4对应风险二级，0.4
‑
0.6对应风险三级，0.6
‑
0.8对应风险四级，0.8
‑
1.0对应风险五级。
[0046]
本发明直接基于sar影像进行洪涝风险区域提取，利用c波段后向散射系数(σ0)计算洪涝事件前后的两时相sar影像的水体概率，再利用两时相影像的水体概率增量表达洪涝风险区域。其中，水体概率为像元属于水体的概率。该方法是一种简单快捷的洪涝监测方法，可以快速自动提取降雨事件后洪涝淹没新增区域(洪涝风险区域)。与现有技术相比，本发明具有如下特点和有益效果：
[0047]
1、本发明具有自动、快捷、高效提取洪涝风险区域的特点，不需要手动选择设置阈
值，利用两时相后向散射系数的混合高斯分布特点，基于降采样数据的迭代计算实现效率提高，以及利用非监督分类、em算法构建水体概率自动化计算方法，并制定基于水体概率的洪涝风险区域提取规则，实现新增洪涝风险区域的自动提取。
[0048]
2、本发明使用水体概率增量表达事件前后每个像元属于水体的概率的变化值，相比以前的二值化淹没区域变化，更能表达复杂水体环境带来的陆地表面水体淹没状况的不确定性以及淹没程度的不确定性，另外，以概率增量的形式表达，洪涝风险区域的可视化效果更强，蕴含内容更丰富。
附图说明
[0049]
图1为本发明提出的一种基于sar影像的洪涝风险区域自动提取方法计算流程图。
[0050]
图2为本发明计算的两时相后向散射系数集合的频率直方图。
[0051]
图3为经本发明处理后两时相sar影像分别对应的水体概率值，其中，左侧和右侧分别显示2020年6月26日和2020年7月8的水体概率值。
[0052]
图4为本发明基于降雨事件前后两时相水体概率计算的洪涝风险区域的空间分布图。
具体实施方式
[0053]
下面结合附图及实施例对本发明做进一步说明。应当注意，这里描述的实施例只用于说明解释本发明，并不用于限制本发明。
[0054]
本发明利用sentinel
‑
1sar数据进行洪涝风险区域提取，两时相影像的观测时间分别为分别是2020年6月26日和2020年7月8日，包括以下步骤：
[0055]
步骤(1)，获取鄱阳湖局部区域内降雨事件发生前后的两时相sentinel
‑
1sar影像，并进行预处理，得10m*10m空间分辨率影像图，将研究区内所有像元组成后向散射系数(σ0)集合σ
10*10
；
[0056]
步骤(2)，对集合σ
10*10
绘制频率直方图，假设分布有水体的两时相sar影像的σ0集合符合混合高斯分布，并对σ0集合中的离群值过滤；
[0057]
步骤(3)，对步骤(2)过滤后对应的两时相sar影像进行5*5降采样至50m*50m，生成σ0集合σ
50*50
，对σ0集合σ
50*50
进行k均值聚类，计算本次实验数据的水体和非水体先验概率；
[0058]
步骤(4)，利用最大期望(em)算法和步骤(3)的水体先验概率对降采样后的集合σ
50*50
进行迭代计算分布参数，以提高迭代速度，然后对每个10m*10m像元进行水体概率计算；
[0059]
步骤(5)，以水体概率增量表达事件前后每个像元属于水体的概率的变化值，确定洪涝风险区与非风险区的水体概率增量边界值，构建基于两时相水体概率的洪涝风险区域提取方法。
[0060]
以下参照图1对本发明进一步详细说明：
[0061]
步骤(1)中，获取区域内洪涝事件发生前后的两时相影像，分别是2020年6月26日和2020年7月8日，预处理后得10m*10m空间分辨率影像，将研究区内所有像元组成后向散射系数(σ0)集合σ
10*10
，具体包括以下子步骤：
[0062]
步骤(11)，获取洪涝事件前后覆盖研究区的两时相sentinel
‑
1干涉宽幅(iw)的地距(grd)产品；
[0063]
步骤(12)，利用sarscape软件对两时相sar影像进行配准、多时相滤波、地理编码、辐射定标、空间裁剪，获取研究区域内垂直
‑
垂直(vv)极化σ0影像(单位：db)，此时，空间参考为wgs84，空间分辨率10m；
[0064]
步骤(13)，将两时相影像所有像元组成σ0集合σ
10*10
。
[0065]
步骤(2)中，对集合σ
10*10
绘制频率直方图，假设分布有水体的两时相sar影像的σ0集合符合混合高斯分布，并对离群值过滤，具体包括以下子步骤：
[0066]
步骤(21)，对σ0集合绘制频率直方图，观察其分布，假设其概率密度曲线符合混合高斯分布；
[0067]
步骤(22)，按数值从小到大排列，认为后0.1％的样本为离群值，将其过滤，对离群值过滤后的频率直方图如图2所示，形成了具有两个高斯分布融合而得的混合高斯分布；
[0068]
步骤(3)中具体包括以下子步骤：
[0069]
步骤(31)，对步骤(2)过滤离群样本后对应的两时相sar影像，也即在去除0.1％离群样本后剩余的两时相sar影像进行降采样至50m*50m，生成σ0集合σ
50*50
。
[0070]
步骤(32)，对σ0集合σ
50*50
进行k均值聚类，由于所有象元σ0值的频率分布直方图呈现混合高斯分布，水体表面较为平滑，具有镜面反射特性，会反射大部分微波信号，因此，水体相对其他地物具有更低的后向散射值，在sar影像上通常为暗黑色。所以，根据地物特点和混合高斯分布可将地物分为水体类和非水体类，σ0值较低的为水体类，σ0值较高的为非水体类；
[0071]
步骤(33)，计算σ0值较低的水体像元集合的像元个数占所有像元的个数比值，该比值为水体先验概率初值，即p0(w)。
[0072]
步骤(4)中，利用最大期望(em)算法和步骤(3)的水体先验概率对降采样后的集合σ
50*50
进行迭代计算分布参数，以提高迭代速度，然后进行每个10m*10m像元的水体概率计算，具体包括以下子步骤：
[0073]
步骤(41)，利用混合高斯分布构建σ0边际分布的概率密度函数；
[0074]
步骤(42)，利用基于最大似然估计的最大期望(em)算法对分布曲线的均值、标准差参数迭代计算，对p(w)进行迭代优化。
[0075]
步骤(43)，将步骤(42)中优化的p(w)、p(σ0|w)和代入水体概率计算公式；
[0076]
步骤(44)，将两时相影像的σ0代入公式，分别计算两时相影像的水体概率(如图3)，其中左侧和右侧分别为该区域2020年6月26日和2020年7月8的水体概率值，概率值阈值为0
‑
1，越趋近于1，该像元属于水体的概率越大，反之，越小。
[0077]
步骤(5)中，将洪涝事件发生前后的两时相影像的水体概率作差，计算洪涝风险区和非风险区的水体概率增量边界值，并以水体概率增量的形式表达结果，具体包括以下子步骤：
[0078]
步骤(51)，求取洪涝事件后前两时相影像的水体概率之差；
[0079]
步骤(52)，0为两时相影像的水体概率之差的关键阈值，利用水体概率进行洪涝风险区提取，并以概率增量的形式表达(如图4)，理论上两时相水体概率之差的阈值为0
‑
1，越
接近于1表明，该像元由陆地变为水体的概率越大。
[0080]
步骤(53)，对y进一步进行阈值分割，设定0.2、0.4、0.6和0.8为阈值，0
‑
0.2对应风险一级，0.2
‑
0.4对应风险二级，0.4
‑
0.6对应风险三级，0.6
‑
0.8对应风险四级，0.8
‑
1.0对应风险五级。