基于局部自适应阈值的GF-3 SAR影像洪涝水体提取

杨雨晴，王冬梅，徐 佳

（1.河海大学 地球科学与工程学院，江苏 南京 210098；2.江苏省水利科学研究院，江苏 南京 210017）

洪涝灾害具有持续时间短、危害性大等特征，一直危害着人类社会的发展。因此，及时防御洪涝灾害，减少灾害损失，保障人类安全，显得至关重要[1]。基于卫星遥感技术，人们可实现快速、准确的洪涝淹没范围提取，直观地显示其空间分布及其动态变化和发展规律，从而使得其在洪涝灾害监测中发挥着越来越重要的作用[2]。经过近几十年的探索和实践，基于遥感技术的洪灾应用，逐步形成了贯穿于灾前、灾中和灾后全过程的洪涝灾害监测。水体信息提取是洪水监测的出发点。传统的监测方法主要通过采样调查，耗时耗力，花费巨大，遥感技术的发展和应用则提供了一种新的选择。光学传感器虽然有着较好的洪水监测潜力，但时间分辨率和空间分辨率往往不能兼顾，且易受云层和天气影响，而合成孔径雷达（SAR）受多云与阴雨天气影响小，可实现全天候、全天时的大面积观测[3]，已成为洪涝灾害监测的重要技术手段。一直以来，我国在民用SAR遥感卫星方面与其他发达国家存在较大差距。2016-08-10成功发射的高分三号卫星是我国首颗C频段多极化高分辨率SAR卫星，其具有高分辨率、大幅宽、多成像模式等特点，缓解了我国SAR数据的短缺问题，为洪涝灾害研究提供了重要的数据源。

对于洪涝灾害遥感监测而言，从遥感影像中快速准确地提取洪涝水体信息至关重要。基于SAR影像的水体信息提取方法[4]主要有阈值分割法、机器学习法、曲线演化法等。由于洪涝灾害监测属于应急监测，更侧重于较高的提取速度，因此阈值分割法在SAR影像洪涝水体提取中应用最为广泛，如李景刚等[5]采用改进的OTSU法进行水体提取研究；郭欣等[6]采用双峰法和OTSU法对Sentinel-1A SAR影像进行洪涝淹没监测等。然而在大场景SAR影像中水体目标相对于整幅影像占比一般较小，全局阈值法往往提取效果不佳，特别是当研究区内存在地形起伏的山区时，山体阴影亦会对水体提取造成干扰。因此，不少学者提出了改进算法，如安成锦等[7]利用多阈值分割方法进行了SAR图像水域分割；Junhua等[8]结合纹理信息，采用非监督K-means聚类获取初始水体范围；庞科臣等[9]利用改进的OTSU法对DEM数据进行分割来去除误提为水体的山体阴影。本文在以上研究的基础上，针对高分三号SAR影像特点，提出了一种DEM辅助下基于局部自适应阈值的SAR影像洪涝水体提取方法，利用GF-3精细条带2成像模式（FSII）SAR影像进行洪涝水体提取实验，并实现了对湖南省东北区域洪水淹没范围的快速提取。

1 研究区和数据介绍

实验区位于湖南省东北部，如图1b所示。该实验区位于益阳市和长沙市交接处，其中，宁乡市所占面积最大。宁乡市，原名宁乡县，2017年撤县建市后，成为湖南省长沙市管辖下的县级市。该实验区介于112°07′E～112°44′E， 28°03′N～28°31′N 之间，属中亚热带向北亚热带过渡的大陆性季风湿润气候，四季分明，温度适宜，雨水充足，寒冷期短，炎热期长。每日平均气温16.8℃，每年平均降水量1 358.3 mm。内部有沩水河及其他支流穿过，有较多山地分布，山间有较小的湖泊和水洼。境内地貌以山地、丘岗、平原为主，平均海拔高度为153 m。

本文实验数据采用一景GF-3 L1影像、SRTM DEM数据和两景Sentinel-2A数据。高分三号数据的获取时间为2017-07-09，成像模式为精细条带2模式，极化方式为HH和HV双极化，距离向分辨率为1.1 m，方位向分辨率是5.8 m，入射角为27.8°，如图1c所示，该RGB图像由R：HV极化；G：HH极化；B：HV极化组成。Sentinel-2A的获取时间为2017-05和2017-07，产品级别为Level-1C，是经正射校正和亚像元级几何精校正后的大气表观反射率产品，其中5月数据为灾前数据，7月为验证数据。

图1 研究区示意图

2 研究方法

为解决大场景SAR影像中水体目标占比较小问题，本文首先将原始影像划分为K×K个局部区域，然后进一步将每个局部区域分为N×N个子块，通过统计子块中初始水体和非水体的比例来筛选子块，从而提高局部区域中水体占比。同时，考虑到山体阴影会造成水体误提和辐射失真，在阈值分割前，利用DEM和SAR轨道参数信息生成阴影掩膜文件从而去除SAR图像的辐射失真区域。本文提出的适用于高分三号SAR影像的洪涝水体提取算法流程如图2所示，主要包含5个步骤：①影像预处理：对GF-3数据进行辐射定标、多视、滤波、地理编码等处理；②山体阴影掩膜文件制作：将DEM和SAR轨道参数信息相结合，生成阴影掩膜文件；③初始水体分布获取：利用山体阴影掩膜文件去除SAR图像的辐射失真区域，然后将K-means聚类方法应用于阴影掩膜后的SAR图像，得到初始的水体分布图；④洪涝水体自适应阈值提取：对影像分区分块处理，利用KI方法对优化后的局部区域进行阈值分割，提取最终的洪涝灾区水体范围。⑤洪涝灾害监测：以Sentinel-2A数据为参考，提取灾前水体分布图，并将本文方法提取的灾后水体分布图与其进行差分处理，得到洪水分布图。

图2 水体信息提取流程

2.1 山体阴影检测

由于SAR成像系统为侧视成像，SAR图像存在固有的几何失真，在地形起伏较大的山区更容易产生。受山体的影响，背向坡将无法获得雷达波束，在影像上就形成阴影，阴影区的大小取决于俯角和斜坡坡度两个参数。雷达阴影在图像相应位置呈现暗色，即使图像经过正射校正，雷达阴影区的信息也无法得到补偿，其后向散射系数和水体的后向散射系数仍然接近，使得图像分割无法有效区别阴影和水体信息。因此，为了消除山体阴影的影响，本文首先采用DEM数据模拟SAR图像，然后根据阴影区的低亮度特性，从模拟SAR图像中检测山体阴影并将其掩膜。

由雷达成像原理可知，DEM模拟SAR图最主要的影响因素是本地入射角的大小，其计算公式[10]为：

式中，θ为雷达入射角；θA、θC为分辨单元在距离向和方位向的坡度角。

在模拟雷达图像时，首先根据分辨率单元的大小将地面分成格网，对地面不同部分的雷达波束入射角进行计算。由于雷达图像为斜距成像，需考虑地距和斜距的关系，即：

式中，R为斜距；G为地距；θ为雷达波束入射角。

根据地物的散射特性和成像原理，按雷达方程和灰度方程，计算每一块地面单元的散射回波强度[11]。雷达方程式如式（3）所示：

式中，pr为接收功率；pt为发射功率；G为天线增益；λ为波长；R为天线和目标直接的距离；σ0为后向散射系数；ΔA为分辨单元。

2.2 KI阈值分割算法原理

陆地水域以镜面散射为主，比植被、城镇等非水体表面的后向散射弱，因此可通过阈值分割方法提取水体信息，当图像的后向散射强度x小于阈值时定义为水体，大于阈值定义为非水体。本文采用KI算法进行图像分割，该方法由Kittler和Illingworth于1986年提出，以图像的概率分布为特征，基于最小错分率和贝叶斯理论来进行阈值选择[12]。

在代价函数的辅助下，当满足总体分类代价最小时即可获得具有最小错分率的结果。

在KI算法中，假设直方图为h(x)，则通过以下目标函数寻找最优阈值：

(2) 专用耐磨轧制钢板NM360系列。NM360是在Q345基础上通过调整合金元素后生产的新钢种，主要应用于矿山机械、煤矿机械、环保机械、工程机械等，也常用作为屈服强度不小于700 MPa高强度结构钢使用。主要是为需要耐磨的场合或部位提供保护。该衬板的屈服强度大于800 MPa，抗拉强度大于1 000 MPa，布氏硬度为360。其主要化学成分见表1。

其中

这种利用控制最小化分类误差而确定阈值的方法可以通过最小化函数而获得[13]，即

式中，T*为最优阈值，即分割阈值。

2.3 基于局部自适应阈值的水体信息提取

考虑到大场景SAR影像中水体目标占比较小，全局阈值法难以获得满意的提取效果。本文通过对原始影像分区分块，在局部区域中利用KI算法计算自适应的局部阈值，实现影像的水体和非水体精确分割。具体步骤如下：

1）提取局部区域。将影像分为K×K个局部区域，K值大小可根据影像大小适当选取，本文K取6。

2）获取新的局部区域。当水体占影像比例较小或较大时，KI方法计算的阈值都不准确，因此为了提高阈值计算的准确性，进一步将每个局部区域分为N×N个子块，根据初始水体分布图对子块进行筛选，当该子块满足水体占图像一定比例（0.5±0.15）时，将这些子块中的所有像元作为新的局部区域。

3）计算局部区域阈值。根据上一步中提取的所有像元值进行判断：当局部区域内几乎完全为水体或非水体时，以全局阈值作为该窗口的阈值；当局部区域有水体与非水体时，将提取的所有窗口的像元值（即新的局部区域）通过KI算法计算出水体和非水体的分割阈值。

3 实验结果分析

3.1 数据预处理

本文使用PIE软件进行GF-3 SAR数据预处理工作，具体包括影像数据导入、辐射定标、复数据转换、多视处理、滤波、地理编码、转dB影像、研究区裁剪等。其中多视处理选择距离向和方位向的视数为5和2；滤波过程中选择7×7窗口的EnLee滤波器进行滤波处理。最终获得空间分辨率为15 m，WGS-84坐标系下的后向散射系数影像，如图3所示。

图3 GF-3 SAR影像后向散射系数影像

3.2 山体阴影掩膜

为了便于后续实验的顺利进行，本研究结合DEM和SAR轨道参数生成了SAR模拟图，如图4所示，可以发现在SAR模拟图中，阴影区域表现为黑色，后向散射系数相对较低，因此通过分割阈值的方法可以提取山体阴影区域，如图5所示，并制作山体阴影掩膜文件，去除山体阴影对水体提取精度的影响。

图4 SAR模拟图

图5 山体阴影提取

3.3 自适应阈值实验结果分析

为了解决全局阈值难以准确提取水体信息的问题，本研究以HV极化为例进行自适应阈值实验研究。首先采用非监督分类K-means聚类方法对山体阴影掩膜后的影像进行分类（K=15），并选择聚类1和聚类2组合为水体类，其他的为非水体类，得到初始的水体范围，如图6所示。从图6可以发现该方法得到的水体分布图存在大量噪声，有部分裸地、土壤等地物的影响，但是由于山体阴影掩膜的加入，山体阴影的影响被有效消除。

图6 初始水体分布图

接下来将山体阴影掩膜后的影像和初始水体分类图分为N×N个局部区域（N的值适影像质量和大小而定，此实验N=6），在初始的水体范围基础上，用10×10的窗口对每个局部区域进行计算，统计像元值满足水体占图像一定比例（0.5±0.15）的所有窗口，将这些窗口中的所有像元作为新的局部区域的水体和非水体，并用KI法进行阈值分割。若局部区域内完全为水体或非水体时，以KI全局阈值作为该窗口的阈值；若局部区域有水体与非水体时，将提取的所有窗口的像元值（即新的局部区域）通过KI算法得到水体和非水体的分割阈值。然后对原始影像进行阈值分割，得到最终的水体分布图，如图7所示。图8为KI方法对山体阴影掩膜后的影像进行全局阈值的结果，可以发现全局统一阈值结果存在部分噪声，且有裸地等地物的误提和少量河流出现断层现象。

图8 全局阈值分割结果

从图7、8可以看出，与全局阈值分割结果相比，自适应阈值分割结果噪声有所降低。区域1为沩水下游区域，如图9a和10a所示。使用全局统一阈值提取水体时，阈值偏小，河流呈现断断续续状态，如图9a红色框中所示，出现了漏提现象。而自适应阈值的水体提取过程中，该现象有所改善，阈值相对较大；区域2为沩水中下游河流，如图9b和10b，全局统一阈值法在水体提取过程中，阈值偏大，由于裸地等的后向散射系数与水体的后向散射系数相近，则会被分为同一类，误提现象明显，如图9b红色框中所示。而自适应阈值法可以有效降低其他易混淆地物的干扰；区域3为七家湾区域如图9c和10c，该区域水体内部被两条大坝所分隔，从全局阈值分割的结果图可以看出，左边的大坝出现断断续续现象，有部分被错分为水体，右边的大坝缺失比较严重，大部分被错分为水体，而在自适应阈值水体提取中，左边的大坝则被完整的分割为非水体部分，右边的大坝的分割结果也有所改善。

图7 自适应阈值分割结果

图9 全局阈值水体提取局部放大结果

图10 自适应阈值水体提取局部放大结果

3.4 水体提取精度验证

本文以灾后7月该地区的Sentinel-2数据为参照数据，将GF-3 SAR数据水体提取结果与Sentinel-2数据的水体提取结果进行比较，采用查全率、查准率和虚警率3个指标对本文提出的自适应阈值水体提取方法和全局单一阈值水体提取方法进行优劣性评价，查全率越高，提取结果可靠性越高，查准率越高，监测准确性越高。随机选取局部区域进行精度评定，统计结果如表1所示。

表1 水体提取精度对比

从表1可以看出，在查全率相差不到4%的情况下，从监测准确性出发，查准率从82.49%提升到90.46%，上升了7.97%。综合考虑查全率和查准率的指数F-Measure，局部自适应阈值方法要优于全局统一阈值方法。

3.5 洪涝灾害监测实验结果分析

该区域在2017-06-10～2017-07-03时间段发生连续降雨现象。本研究对2017-05-18的Sentinel-2A影像和2017-07-09的影像分别进行了水体提取，获得灾前和灾后水体，再将灾前和灾后水体叠加，获得最终的洪水分布图。比较该区域的灾前灾后的水体叠加图，可以发现团头湖附近有明显的水体变化情况，如图11所示。

图11 团头湖洪水分布图

4 结 语

本研究基于高分三号数据，提出了一种有效的地形数据支持下的自适应阈值洪水水体提取方法。本文以湖南省东北部地区为研究对象展开实验，主要结论如下：

1）该研究区水体面积占整幅影像的比例较小，KI法计算的阈值偏大。利用自适应阈值水体提取的方法可以有效获取水体信息，去除部分裸地、土壤等地物的影像，并保留了水体的细节信息。

2）在SAR影像中阴影信息必不可少，且会对水体的提取造成一定的干扰。本文加入了地形数据，去除了辐射失真区域，有助于获得更准确的水体-非水体信息。

3）从该实验区的洪涝灾害监测的实验结果来看，该方法对快速、准确的洪涝灾害监测具有很好的适用性。

本文方法也有不足之处， 由于KI法计算的全局统一阈值偏大，因此保留了河流边缘的细节信息，导致噪声和相似地物的干扰也会增强，而本文的方法虽然大大减少了噪声和相似地物的影响，但是有些地方的水体细节保留情况不如KI法计算的全局统一阈值方法，使得查全率低于KI法计算的全局统一阈值的查全率。