基于Sentinel-1数据的典型地物特征分析与洪水区面积提取

李 颖， 陈怀亮， 叶昊天， 周兆基

(1.中国气象局·河南省农业气象保障与应用技术重点开放实验室，郑州 450003；2.河南省气象科学研究所，郑州 450003； 3.哈尔滨市气象局，哈尔滨 150028； 4.河南省气象局，郑州 450003)

引 言

近年来，全球气候变化加剧[1]，极端天气现象频繁出现。2021年7月17日至22日,河南省出现了历史罕见的极端暴雨过程[2]。由于强降水引发的洪涝灾害发生频率较高、影响范围较大，因此在灾后快速获取大面积洪水分布信息对于灾害评估、救援等具有重要意义。遥感影像具有宏观、客观、及时获取地表信息的优势，可以克服灾情现场调查时交通不便和淹没区域难以抵达的困难，节约人力成本和时间成本，在洪水面积提取方面具有重要的应用价值。

晴空条件下中高分辨率光学遥感影像在水体信息提取方面具有较高精度，已得到广泛应用[3-4]。由于光学影像难以穿透云层，而洪水发生前后通常为云雨雾天气，因此使用光学遥感影像难以及时获取地表信息。合成孔径雷达(SAR)具有穿透云层的能力，并可通过地物在不同极化方式下的后向反射特性识别地表类型，在灾害性天气发生过程中可以全天时、全天候持续监测地表变化情况[5]。现有研究中，将SAR影像用于水体信息提取的技术已较为成熟，最广泛应用的方法主要有基于极化灰度图、干涉相干图、纹理图等的阈值法[6-8]和自动聚类算法[9-10]。相关研究包括Marti-Cardona等[11]通过入射角分析识别湿地洪水区和洪水造成的开放水体；孙亚勇等[12]通过山体阴影剔除，提高了山区水体提取的精度；谷鑫志等[13]在阈值分割的基础上，综合利用了GF-3影像的多极化信息与空间上下文信息，利用最大后验概率准则输出水体分布图等；眭海刚等[14]利用深度特征和语义信息实现了灾前光学影像和灾后SAR影像的高精度配准,进而实现洪水变化监测和灾损信息提取；严夏青等[15]基于Sentinel-1数据提取长江干流水域面积并开展动态变化监测。

总体来看，在基于SAR影像的水体识别和洪水面积提取方面，国内研究成果相对较少。一方面，对包含水体在内的典型地物进行后向散射特征分析时，多数研究通常仅针对不同极化方式进行分析，对局地入射角对典型地物后向散射特性影响的研究较少；另一方面，在洪水发生后往往仅对新增开放水体进行识别，没有开展新增湿地的识别，故提取洪水区面积偏小。针对上述问题，本研究统计分析了不同典型地物的后向散射系数随局地入射角的变化情况，并在洪水发生后提取新增水体的基础上，进行了新增湿地与易混淆地物的极化分解分析，进而提取了新增湿地，综合新增开放水体与新增湿地得到更为完整的洪水区面积。

1 研究区与数据源

1.1 研究区

美国2018年的Florence飓风为30年来威胁南北卡罗来纳州的最强级别飓风。该飓风于2018年9月13日登陆美国东海岸，随后给南北卡罗来纳州带来了强降水和灾难性洪水，其洪水至9月26日仍未完全消退。本文研究区为覆盖南北卡罗来纳州部分地区的一景Sentinel-1影像，包含Catawba河和Congaree河两条由Florence飓风造成的洪水泛滥的河流的部分区域。

1.2 数据源及预处理

1.2.1 Sentinel-1数据

Sentinel-1载荷为C波段SAR传感器，工作频率为5.4 GHz，由Sentinel-1A和Sentinel-1B两颗卫星组网，单颗星重访周期为12 d，双星组网的重访周期为6 d。Sentinel-1的成像模式包括条带、干涉宽幅、超宽幅和波浪。本文共使用了6幅Sentinel-1A和Sentinel-1B干涉宽幅模式的level-1级地距影像(Ground Range Detected，GRD)产品，该产品采用WGS84椭球将侧视成像的斜距图像校正为地距影像，使影像的方位向和距离向的分辨率一致，均为10 m，具有VV和VH两种极化方式。根据数据质量和轨道情况(选择同轨道数据)，影像成像时间为2018年7月26日、8月7日、8月19日、8月31日、9月12日和9月18日，包含洪水前的5景影像和泛滥期的1景影像，以支持研究区Florence飓风造成的洪水淹没遥感监测。

1.2.2 雷达数据预处理

使用SNAP软件和ENVI软件对Sentinel-1 level-1级地距影像进行预处理，主要步骤包括热噪声消除、多视处理、辐射定标、斑噪去除、地形校正和直方图拉伸。其中，多视处理的等效视数设为4.9，辐射定标得到Sigma后向散射系数，斑噪去除选择Gamma-Map滤波器(滤波窗口设为5×5)，地形校正采用SRTM DEM数据(空间分辨率为90 m)。详细的处理方法见参考文献[16-17]。预处理后斑点噪声得到较好抑制，水域边界清晰光滑，9月18日的影像较洪水前的2景影像可见明显的泛滥区(图1)。

图1 预处理后研究区局部地区VV极化、VH极化、VV极化RGB合成图

1.2.3 NLCD2016数据及预处理

本文使用美国地质调查局(USGS)发布的2016年国家土地覆被数据库(NLCD)分类数据(NLCD2016)作为洪水前的地表真实覆被数据，用于SAR图像中典型地物特征分析与解译。裁剪出覆盖研究区范围的NLCD2016数据，根据研究区实际情况对地物类别进行合并，归纳为水体、裸土、建成区、农作物、灌木、森林和湿地等7种典型地物(图2)。

图2 研究区土地利用类型分类图

2 典型地物特征分析

2.1 后向散射系数均值分析

利用NLCD2016数据对洪水前的5景影像进行分类，计算不同典型地物在两种极化方式下后向散射系数的均值(表1)。由表1可看出，各种地物类型VV极化的后向散射系数均高于VH极化的后向散射系数，同时VV极化后向散射系数较高的地物VH极化后向散射系数也较高，反之亦然。两种极化方式下，水体后向散射系数的均值均最低，其次较低的是裸土的，最高的则是湿地的。在VV极化下，建成区、农作物、灌木、森林等均值的可分性较差；在VH极化下，森林后向散射系数的均值则明显高于另外三种地物的均值。

表1 典型地物的后向散射系数均值 dB

2.2 局地入射角分析

在不考虑角度效应的情况下，使用VV和VH两种极化方式，水体、裸土、森林、湿地的可分性相对较好。但考虑地形起伏后，雷达波的局地入射角(Local Incidence Angle, LIA)亦会在一定范围内变化。LIA是影响地物散射特性的一个重要因素。本文统计不同典型地物的后向散射系数随LIA的变化情况，结果显示，LIA越小，地物的回波强度越大；当LIA在0°～90°均匀变化时，所有地物类型在VV和VH极化方式下均无法仅使用后向散射系数进行有效区分(图3)。

图3 典型地物的后向散射系数随局地入射角的变化情况

实际上，大部分星载雷达的入射角设计在40°左右，因此雷达波的LIA均匀变化在0°～90°的可能性较小。本文统计了研究区6景雷达影像中LIA的真实分布情况，统计结果显示，各景影像中90%左右像元的雷达波LIA分布在30°～45°(图4)。分析LIA在30°～45°时典型地物的后向散射系数特征可见，水体在两种极化方式下较其他地物均有较好的可分性，特别是在VV极化方式下，使用后向散射系数可有效识别开放水体。

图4 局地入射角在30°～45°时典型地物

2.3 纹理特征分析

现有基于雷达影像的地物分类、识别研究中，纹理特征使用较为广泛。本文亦使用灰度共生矩阵(GLCM)分析研究区各类典型地物的多种纹理特征。GLCM是一种通过图像灰度空间的相关性描述纹理的方法。其定义为：在图像内，统计所有从灰度级为i的像元点，沿移动方向移动一定距离后，到达灰度级为j的像元的概率为p(i,j)时所形成的方阵。不同地物类型的7种GLCM纹理特征，包括方差(Variance)、均质性(Homogeneity)、对比度(Contrast)、熵(Entropy)、角二阶矩(ASM)和相关性(Correlation)，计算公式如下：

(1)

(2)

(3)

(4)

(5)

(6)

其中，N为图像的量化等级，μ为GLCM的均值，σ为GLCM的标准差，其计算公式如下：

(7)

(8)

本文计算GLCM纹理特征时选择3×3的窗口，移动方向设为0°、45°、90°和135°，图像量化等级设为64级，步长设为1。分别统计7种典型地物在VV极化方式下的6种纹理特征(表2)，结果显示：建成区在方差和对比度上明显大于其他地物类型的；水体则在6种纹理特征上均区别于其他地物类型的，特别是方差与对比度明显小于其他地物类型的；在均质性、熵和角二阶矩上，裸土和建成区的均值接近，而农作物、灌木、森林和湿地的均值接近。由以上分析可知，纹理特征可有效应用于水体识别，并可帮助区分植被与非植被。

表2 典型地物的6种纹理特征

3 洪水区域提取方法

3.1 ISODATA算法

ISODATA算法属于非监督分类方法，是在K-均值算法的基础上改进的迭代自组织聚类方法，该方法在开放水体的提取中已得到广泛应用[18]。通过上文分析，使用VV和VH两种极化方式下的后向散射系数与方差和对比度两种纹理特征，可有效识别开放水体和其他地物类型，故选择上述4种特征值进行ISODATA迭代分类。分类后进行类别合并，将水体作为一类，其他地物类型合并为一类，之后运用数学形态学算子，通过聚类处理消除开放水体中的斑点。

3.2 极化分解

发生洪水后，一方面包括裸土、道路、农作物、牧草等在内的低矮地物可能被淹没形成开放水体，另一方面森林和灌木的底部可能被淹没，从而形成湿地。利用Freeman极化目标分解理论分析洪水发生后由森林和灌木变化形成的湿地区域的后向散射特征变化。Freeman分解由Freeman和Durden于1998年提出[19-20]，该方法在假定散射体满足反射对称性的条件下，对极化协方差矩阵建立包括适度粗糙表面的Bragg散射项、随机方向偶极子组成的体散射项和由正交平面构成的偶次散射项等3个分量的基本散射机制模型。Bragg散射的典型地物包括水体、裸土、硬化地面等，体散射的典型地物包括森林、灌木等，偶次散射的典型地物则包括复杂建筑物、建筑与地面、冠层与地面、冠层与水面等能够形成角散射的地物组合。其中，冠层与水面形成角散射的情况主要发生在湿地。

进行Freeman分解时，森林和灌木的主要散射分量为体散射，另外两种分量贡献较小，在非定量研究中较之体散射可忽略。湿地则存在来自森林、灌木冠层的体散射和冠层与水面之间的偶次散射两种主要散射分量，且其偶次散射项的回波信号相比体散射的后向散射系数较高。因此，针对洪水发生前后的森林和灌木，消除回波信号中的体散射项后，若仍存在较高的回波信号，则说明发生了偶次散射，即被洪水淹没。Freeman分解中体散射的极化协方差矩阵形式为：

(9)

式中，fv表示体散射分量的权重。

在Freeman分解理论中，体散射的同极化后向散射系数是交叉极化的3倍。因偶次散射项理论上不存在交叉极化，因此应用Freeman极化目标分解理论，计算极化分解系数

σ3v=SVV-3SVH[21-22]

(10)

其中，SVV为VV同极化后向散射能量，SVH为VH交叉极化后向散射能量。通过计算σ3v，体散射项可以得到较大程度的消除，而偶次散射项不受影响，即森林与灌木的极化分解系数σ3v很低，而湿地的σ3v较森林与灌木明显较高。表3为洪水发生前后森林、灌木和湿地的极化分解系数。由表3可看出，洪水发生前5幅影像中森林的平均σ3v最低，其次是灌木的，湿地的平均σ3v明显高于前两者的。洪水发生后，原森林区域和灌木区域的平均σ3v明显升高，显示研究区域中确有大面积的灌木和森林因洪水淹没变为湿地。在上述分析的基础上，通过判读洪水前森林、灌木、湿地的极化分解系数的变化特征，结合VV后向散射系数设置变化检测阈值，可提取洪水后形成的湿地。

表3 洪水发生前后典型地物的极化分解系数

4 结果与分析

4.1 开放水体提取

利用ISODATA算法和VV、VH两种极化方式下的后向散射系数与方差、对比度两种纹理特征，分别提取洪水发生前雷达影像中的开放水体和洪水发生后雷达影像中的开放水体，计算其差值并做去噪声点等后处理，即得到洪水发生后新增开放水体面积(共计449938个像元，新增水体面积45 km2)，局部提取效果见图5，图例与图2相同。洪水淹没区域主要分布在Catawba河沿岸，Congaree河沿岸也有部分区域被淹没。对照NLCD2016土地覆被数据，淹没区域的土地利用类型主要是农作物、灌木和湿地。

图5 洪水后新增开放水体提取结果

4.2 新增湿地提取

在分析洪水发生前后森林、灌木、湿地的极化分解系数变化特征的基础上，排除洪水发生前的湿地面积，利用阈值判别法和VV极化方式下的后向散射系数提取新增湿地。因新增湿地区域出现在洪水前的森林和灌木土地覆盖类型上，通过典型地物雷达回波特征分析，这两类土地覆盖类型在变化为湿地后，VV极化后向散射系数将明显提升。结合图3与图4中典型地物的后向散射系数随局地入射角的变化和分布情况，在极化分解原理的基础上，经人机交互实验，提取洪水发生后森林和灌木VV极化后向散射系数值增大3 dB以上的区域为新增湿地面积，做去噪声点等后处理，共计得到2124106个像元，面积约212.4 km2，局部提取效果见图6，图例与图2相同。新增湿地区域主要分布在现有湿地的周边。

图6 洪水后新增湿地提取结果

5 结论与讨论

本研究将洪水发生后淹没区分为两种情况：一种是发生洪水后，裸土、道路、农作物、灌木、湿地等在内的低矮地物被淹没形成的开放水体;另一种是发生洪水后，森林和灌木下部被淹没形成新的湿地。通过典型地物特征分析，开放水体在VV和VH两种极化方式下较其他地物具有良好的可分性，其方差与对比度两种纹理特征也明显小于其他地物类型的。利用ISODATA算法进行非监督分类提取水体，结果显示，洪水淹没区域面积约45 km2，主要分布在Catawba河沿岸和Congaree河沿岸，即新增开放水体出现在原有开放水体的邻近位置，分析原因可能是河道水位升高，进而洪水漫灌或因河谷地带地势较低容易积水，洪水汇聚形成开放水体。

通过典型地物特征分析，湿地在VV极化方式下后向散射系数的均值明显高于其他地物类型的均值，但其纹理特征与其他地物类型，特别是与森林和灌木的可分性较差。采用极化分解分析显示，洪水发生后原森林区域和灌木区域的平均极化分解系数明显升高，表明研究区域中确有大面积的灌木和森林因洪水淹没变为湿地。经人机交互实验，确定洪水发生后森林和灌木区域VV极化后向散射系数值增大3 dB以上的区域为新增湿地区域，提取新增湿地面积约212.4 km2。

本研究较之国内现有研究往往仅进行新增开放水体的提取，且较少考虑入射角对不同典型地物的后向散射系数影响的问题，在利用ISODATA算法提取开放水体的基础上，分析了不同典型地物的后向散射系数随局地入射角的变化情况和不同地物类型的纹理特征，进行了洪水发生后新增湿地与灌木、森林的极化分解分析，利用阈值法提取了新增湿地面积，综合新增开放水体与新增湿地得到更为完整的洪水区面积。因本研究采用历史卫星遥感影像研究了国外的一次洪水过程，难以开展实地灾情调查，故利用目视解译分析评价本方法提取洪水区位置的合理性。下一步的研究中，拟利用本方法对国内洪水过程开展洪水区面积的实时提取，并开展实地调查，进一步评价本方法的精度与适用性。