孙 诚, 沈 芳, 唐儒罡

（华东师范大学 河口海岸学国家重点实验室, 上海 200241）

0 引 言

我国长期以来饱受洪涝灾害的频繁侵扰, 严重威胁了人民群众的生命和财产安全. 遥感技术不受空间限制, 且可迅速地获得洪水淹没信息, 已成为洪涝灾害监测和评估的常用手段[1-2]. 星载SAR(synthetic aperture radar, 合成孔径雷达) 工作于微波波段, 具有地表穿透能力, 不受云雾限制, 突破了光学遥感技术受天气影响大的局限性[3], 在洪水监测方面具有独特的优势.

阈值分割法及机器学习方法常被用于SAR 影像水体信息的提取. 阈值分割法的关键在于在图像灰度的分布中确定合适的阈值, 以区分水体和非水体. 贾诗超等[4]研究了Sentinel-1SAR 双极化数据之间水体信息提取的关系, 提出了基于阈值分割的SDWI (sentinel-1 dual-polarized water index) 水体指数法. 阈值分割法的操作逻辑简单, 计算时间短, 但在地物类型复杂的情况下, SAR 影像上灰度值相等的像元并不能完全对应同一地物类型, 难以选择最优阈值, 使得阈值方法易受到图像噪声和强度不均匀性的影响, 因此基于影像单个像元的方法具有局限性[5].

为了更好地获得SAR 影像中的信息, 有学者采用GLCM (gray-level co-occurrence matrix) 提取SAR 影像中的纹理信息, 建立多特征空间, 并利用机器学习模型进行水体信息的提取[6-7]. Lyu 等[8]提出了结合4 种灰度共生矩阵纹理特征与SVM (support vector machines) 分类器, 在区分水体目标区域与其他地物及对灰度共生矩阵提取水体等方面进行了初步探索. 胡德勇等[9]使用Radarsat-1 单波段单极化数据对水体和居民地信息进行了提取, 但由于单极化影像信息有限和分类算法的限制, 存在一定的错分、误提. 此外, 与单极化SAR 影像相比, 全极化SAR 影像的信息更为丰富, 图像分类的性能更高[10]. 也有学者采用RF (random forest)[11]、GBDT (gradient boosting decision tree)[12]等机器学习算法, 但传统的机器学习算法由于内存占用较大, 存在模型训练时间较长的问题[13]. 近年来, 基于集成学习的LightGBM (light gradient boosting machine, 轻量级梯度提升器) 算法, 因其学习能力强及预测精度高的特点, 越来越多地被应用于各类学科领域[14-16], 相比于SVM、RF 及GBDT 等机器学习算法, 其具有更快的训练速度和更低的计算代价, 可以快速地处理大数据量和高特征维度的数据[17]. 因此, LightGBM 算法更适用于对速度与精度有较高要求的洪涝灾害淹水信息应急提取.

本文基于Sentinel-1SAR 影像, 旨在建立一种结合SAR 纹理信息和LightGBM 算法的洪水淹没地区遥感应急监测方法, 并将该方法的水体提取结果与SDWI 水体指数法、SVM、RF 及GBDT 等其他方法进行定量和定性对比, 以测试该方法的提取精度和运行效率, 最后对淮河流域中的典型洪水淹没区域进行淹水范围提取和分析应用实践.

1 研究区与数据

1.1 研究区概况

淮河流域是我国南北区域的自然分界线, 降水时空分布不均, 每年的6—9 月为汛期, 10 月至次年5 月为非汛期, 年平均降水量为920 mm, 汛期年降水量最高达1 600 mm, 流域平均每5 年就会发生一场较为严重的洪涝灾害. 因此, 本文的研究区域选取了淮河流域中游主干道附近地区, 该区域内支流众多, 源短流急, 在汛期易形成洪涝灾害. 研究区域范围见图1, 底图为Sentinel-1 影像.

图1 研究区地理位置Fig. 1 Map of the study area

1.2 数据及其预处理

采用SAR 影像为IW (interferometric wide swath) 模式下Level-1 级别Sentinel-1 的GRD (ground range detected) 影像, 包括VH (垂直发射水平接收)、VV (垂直发射垂直接收) 两种极化模式, 分辨率为10 m, 重访周期最短为6 d. 利用欧空局开发的SNAP 软件, 对该模式下的影像进行了轨道校正、热噪声去除、辐射定标、相干斑滤波、分贝化处理及地理编码等预处理. 轨道校正对从欧空局网站下载的Sentinel-1 原始影像中的原始轨道数据进行校正; 热噪声是SAR 系统自带的背景噪声, 会影响雷达得到后向散射信号的精度, 需要在软件中进行热噪声去除; 辐射定标是将SAR 接收的后向散射信号转化为后向散射系数; 相干斑是SAR 影像分类时的噪声, 在软件中过滤相干斑可提高影像分类的精度; 由于SAR 接受的后向散射信号之间差距不明显, 使用分贝化处理的方式可指数化后向散射系数,便于可视化和影像分类; 最后使用地理编码的方法对SAR 影像进行地理坐标的校正, 完成预处理[18].

研究区域内选取成像时间为2020 年3 月5 日、2020 年3 月17 日、2020 年5 月16 日及2020 年11 月12 日共4 天的Sentinel-1SAR 影像作为数据源. 使用随机抽样的方法在研究区域内裁剪样本区,并利用成像时间相同的总体云量小于20%且水体周围无云干扰的Sentinel-2 影像, 对样本区进行目视解译, 获得样本区的水体范围信息. 最终共得到训练集552 万个, 验证集125 万个.

2 研究方法

本文提出了一种基于SAR 纹理信息和LightGBM 算法的WEGL (water extraction based on GLCM and LightGBM) 方法进行水体提取. WEGL 方法由两个主要步骤构成: 第一是SAR 影像纹理信息的提取 (详见2.1 节); 第二是LightGBM 算法构建, 完成LightGBM 参数优化和纹理特征变量优化 (详见2.2 节). 首先对研究区域的Sentinel-1 影像进行预处理和裁剪, 利用Sentinel-2 光学影像与Sentinel-1 SAR 影像进行配准, 得到样本区的水体分布, 进一步基于GLCM 计算SAR 纹理信息, 构建训练样本数据, 进行LightGBM 算法训练, 得到水体提取结果. 具体步骤实现的流程见图2.

图2 方法实现流程Fig. 2 Flow chart of the method performance

2.1 基于GLCM 的SAR 影像纹理信息提取

GLCM 最早由Haralick 等[19]提出, 是一种描述在预定计算窗口内相邻像元或间隔一定距离像元灰度相关关系的矩阵, 常使用于地物复杂的遥感影像中[20-21]. 本文试验使用的8 个纹理特征及其描述见表1. 表1 各式中:N是GLCM 窗口的行列数;i,j为图像灰度级数;u代表均值;σ代表方差;Pi,j代表从图像灰度为i的像素出发, 灰度为j的像素出现的概率.

表1 纹理特征表达式及描述Tab. 1 Expressions of texture feature and their description

影像的纹理信息与计算灰度共生矩阵选择的方向、步长、影像的灰度级及窗口的大小有关[22-23].经过大量试验, 选择0°、45°、90° 及135° 这4 个方向的GLCM 平均值, 步长大小选择1, 影像灰度级选择64, 窗口大小选择7 × 7, 对表1 中的8 个纹理特征信息进行提取.

2.2 LighGBM 算法构建及优化

LightGBM 是一种基于GBDT 框架的分类模型, 其带有深度限制的Leaf-wise 叶子生长策略及直方图作差加速等创新技术[24], 可支持高效的并行训练, 能够快速处理大数据量和高特征维度的数据[25].

2.2.1 精度评价指标

选取了总体精度、交并比、F1指标及Kappa 系数共4 种精度评价指标, 对试验结果精度进行定量分析.

总体精度 (AO) 常用于衡量图像分类的性能, 表示预测正确的水体像元和陆地像元占总像元数量的比重, 其表达式为

式(1)中:nTP表示水体像元被判定为水体像元的数量;nFP表示陆地像元被判定为水体像元的数量;nTN表示陆地像元被判定为陆地像元的数量;nFN表示水体像元被判定为陆地像元的数量.

交并比 (intersection of union, IoU) 表示预测结果与真实地物的相关度, 其表达式为

F1指标是衡量精准率和召回率的综合指标, 其表达式为

式(3)中: 精准率 (P) 表示真实水体像元占所有预测为水体像元的比例, 表达式为

召回率 (R) 表示真实水体像元占所有实际为水体像元的比例, 表达式为

Kappa 系数 (K) 为衡量分类精度的经典指标, 其表达式为

式(6)中:p0代表水体和陆地像元中正确分类的像元数量,pe的表达式为

式(7)中:a1和a2分别代表水体和陆地的实际像元数;b1和b2分别代表水体和陆地的预测像元数.

2.2.2 参数优化

在进行训练之前, 需要为LightGBM 算法寻找最优参数, 确定最优识别精度下的最优参数组合.在其他参数都保持默认参数的情况下, 使用遍历的方法寻找LightGBM 算法收敛的迭代次数. 经过试验, 200 次迭代之后, LightGBM 算法的水体识别准确率趋于稳定, 不再有明显的提升. 保持LightGBM算法迭代次数为200 次, 进一步测试算法识别准确率随树的最大深度 (max_depth) 变化的情况, 试验结果如图3 所示. 当树的最大深度较小时, LightGBM 算法的准确率曲线随迭代次数收敛的速度变慢,且识别准确率较低. 随着树的最大深度的提高, LightGBM 算法的准确率曲线随迭代次数收敛的速度加快, 且识别准确率明显升高. 从图3 中可看出, 树的最大高度为8 时, 模型的识别准确率表现最好.在确定模型最佳迭代次数和树的最大深度后, 使用网格滚动调参的方式, 构建不同参数的排列组合,代入模型中进行测试, 以得到最优识别精度下的最优参数组合. LightGBM 算法的相关参数和调参的范围及步长如表2 所示.

表2 LightGBM 算法参数及调优Tab. 2 Parameters and optimization of LightGBM algorithm

图3 LightGBM 算法识别准确率Fig. 3 Accuracy of LightGBM algorithm for identification

2.2.3 纹理特征变量优化

通过采用LightGBM 算法的信息增益函数, 表征纹理特征变量的重要性排序, 以代表特征变量对模型训练的贡献程度[14]. 图4 分别显示了VH、VV 两组单极化影像中, 各纹理特征变量对模型训练的重要性排序及对识别水体的贡献度. VH 影像中, 特征变量重要性排序前4 位的特征分别为: 均值、均质性、相异度、相关性; VV 影像中, 特征变量重要性排序前4 位的特征分别为: 对比度、均值、相关性、角二阶矩. 进一步, 将上述两组特征重要性排序前4 位的特征组合成8 个纹理特征并进行模型训练, 得到的LightGBM 算法水体提取精度为98.40%, 与使用16 个纹理特征得到的98.41%精度基本一致. 因此, WEGL 方法剔除了一些特征贡献度较低的纹理特征, 在保持精度基本一致的情况下, 可缩减一半的数据量和内存占用, 大大提高了水体提取效率.

图4 VH、VV 极化影像特征重要性Fig. 4 Importance degree of VH/VV SAR image features

3 结果与分析

3.1 多方法对比及精度评价

Sentinel-1 系列卫星的特点是有多种成像方式, 可实现单极化、双极化等不同的极化方式[26]. 为验证基于双极化SAR 影像的分类性能是否优于单极化影像, 将单极化影像和双极化影像加入对比试验中. 为验证LightGBM 算法的分类性能, 在保持训练集相同的情况下, 与LightGBM 同时进行了SVM、RF 及GBDT 算法的预测试验, 并记录各方法在不同数据特征条件下的预测精度及训练耗时,具体参见表3.

表3 不同算法的对比试验Tab. 3 Comparative tests of different algorithms

从表3 的试验结果看出, 对于每种算法而言, 基于VH + VV 双极化方式的SAR 影像的算法提取精度要优于VH、VV 单极化方式的SAR 影像, 说明双极化SAR 影像中的纹理信息更为丰富, 水体识别的性能更高, 故WEGL 方法选择双极化SAR 作为输入的影像. 在水体提取精度方面, 基于双极化影像的SVM 的水体提取准确率为95.19%, RF 和GBDT 的准确率均高于SVM, 分别为96.23%、96.41%. LightGBM 算法的水体提取准确率最高, 为98.40%. 在训练耗时方面, SVM 的用时最长, 共用时超过4 h, 说明在特征维度较大的情况下, SVM 算法拟合时间较长. RF 和GBDT 算法的训练耗时有进一步优化, 为1 h 左右, 基于决策树思想的算法的训练效率明显提升. LightGBM 算法的训练耗时最短, 仅为3 min, 说明在处理海量数据上LightGBM 算法具有明显优势, 更适用于洪水淹没地区应急提取等场景.

3.2 水体提取与精度

为验证WEGL 方法的合理性及有效性, 在研究区域内选取了河道区、湖泊区及洪水淹没区3 类典型区域进行方法的评价, 并将其与SDWI 水体指数法、SVM、RF 及GBDT 算法进行对比.

图5 区域为河道区, 经目视解译和百度地图验证, 1 处为含水量较高的浅滩, 2 处为河道边的道路.SDWI 水体指数法对河道提取较为完整, 但存在错分现象, 其将1 处的浅滩和2 处的道路识别为水体;对比前者, SVM 算法在1、2 处仍存在错分现象, RF、GBDT 算法在2 处无错分现象, 但将1 处的浅滩识别为水体; WEGL 方法识别出了完整的河道部分, 在1 处和2 处的错分现象明显改善, 表现最好.

图5 河道区水体提取结果Fig. 5 Water body extraction results of river

图6 区域为湖泊区, 3、4、5 处为湖泊水体的边缘. 水体边缘处水深较浅, 与水体内部的像元灰度值存在差异. 基于像元的SDWI 水体指数法无法准确识别边缘水体, 存在较多错分和漏提取现象;SVM、RF、GBDT 算法去除了大量噪声影响产生的错分现象, 但在3、4、5 处水体边缘的漏提取现象仍然存在; WEGL 方法对水体边缘的提取更完整和平滑, 错分和漏提取的现象较少.

图6 湖泊区水体提取结果Fig. 6 Water body extraction results of lake

图7 区域为洪水淹没区, 6、7 处为泥沙含量较高的洪水水体. SDWI 水体指数法提取的水体在6、7 处内部存在噪声现象; SVM 算法在6、7 处改善了水体提取的噪声, 但漏提取现象仍存在; RF、GBDT 算法的错分现象较少, 但仍存在水体内部噪声; WEGL 方法对洪水水体的提取结果表现最好,水体内部噪声和漏提取现象较少.

图7 洪水淹没区水体提取结果Fig. 7 Water extraction results of flooded area

各方法提取水体精度的定量评价结果如表4 所示. WEGL 方法的相关评价指标最好, 与SDWI 水体指数法比较, 总体精度、IoU、F1指标及Kappa 系数分别提高了2.00%、16.89%、0.11 及0.12; 与SVM 比较, 总体精度、IoU、F1指标及Kappa 系数分别提高了0.60%、5.35%、0.04 及0.04; 与RF 比较, 总体精度、IoU、F1指标及Kappa 系数分别提高了0.34%、6.10%、0.04 及0.04; 与GBDT 比较, 总体精度、IoU、F1指标及Kappa 系数分别提高了0.24%、3.30%、0.02 及0.02.

表4 水体提取结果与精度对比Tab. 4 Accuracy comparison of water extraction results

3.3 2020 年洪水灾害遥感监测实例

使用WEGL 方法对2020 年7 月淮河流域洪水受灾较严重的霍邱和颍上段进行洪水淹没监测. 图8分别展示了霍邱和颍上段灾前灾中的水体范围及淹没范围的空间分布. 监测表明, 2020 年7 月27 日淮河流域霍邱和颍上段北岸受灾情况严重, 相比于灾前(2020 年7 月3 日)的监测影像, 水面积由194.92 km2变为343.50 km2, 监测区域内的淹没面积达到148.58 km2. 根据霍邱县人民政府网站(http://www.huoqiu.gov.cn) 和颍上县人民政府网站 (http://www.ahys.gov.cn) 的报道, 2020 年7 月两县受强降水影响严重, 霍邱县农作物受灾面积为542.27 km2, 颍上县湖区受灾面积共17.39 km2, 其中农作物为15.83 km2. 上述受灾统计结果较为碎片化, 无法精确到具体的受灾区域, WEGL 方法的洪水淹没监测结果更为客观和完整, 有助于对受灾情况的整体评估.

图8 霍邱和颍上段灾前灾中水体信息提取情况Fig. 8 Outline drawing of water extraction results in pre-flood and flooding of Huoqiu and Yingshang County

4 结论与展望

针对目前常用的SAR 影像水体提取方法精度不足、运行效率低等问题, 通过GLCM 提取SAR 影像纹理信息, 并进行LightGBM 算法参数优化和纹理特征变量的优化, 建立了WEGL 方法, 大幅提高了水体提取的运行效率, 并在水体提取精度上有进一步优化. 经试验分析, 与SDWI 水体指数法、SVM、RF 及GBDT 算法相比, WEGL 方法提取河道、湖泊和洪水淹没区的精度均具有优势, 在一定程度上抑制了道路、裸地等地物的影响, 提取的水体边缘更加清晰且完整. 除了目标提取精度的优势,WEGL 方法的运行效率也显著提升, 更加适用于洪涝灾害淹没地区的应急监测. 将WEGL 方法应用于淮河流域霍邱和颍上段的洪涝灾害监测, 发现2020 年7 月洪水期间, 水面积由194.92 km2变为343.50 km2, 受灾面积达到148.58 km2, 结果表明WEGL 方法具有时空可扩展性, 可用于不同时期和区域的洪涝灾害监测.

WEGL 方法成功实现了洪涝灾害期间淹水范围的快速监测, 为后续研究中需要大面积水域信息快速提取的场景提供了新的思路. 同时, 可进一步开发集成WEGL 方法的软件系统, 实现淹水信息和洪涝受灾情况的自动获取. 因受限于Sentinel-1 卫星的重访周期, 对重点地区进行全天候监测的难度较大. 在后续研究中可尝试基于WEGL 方法, 加入不同类型的SAR 数据, 提高卫星的观测频率, 进一步提高洪涝灾害淹没范围监测的时效性.