基于Sentinel-1B SAR数据的2018年寿光洪水遥感监测及灾害分析

栾玉洁，郭金运，高永刚，刘 新

(1.山东科技大学 测绘与空间信息学院，山东 青岛266590； 2.福州大学 环境与资源学院，福建 福州 350116)

星载合成孔径雷达(SAR)工作于微波波段，具有一定的地表穿透能力，不受云雾限制，在洪水监测方面具有独特的优势，受到了许多学者的重视。Shu等[1]、汤玲英等[2]分别基于SAR数据进行洪水监测，指出了雷达传感器在洪水监测方面的独特优势，验证了利用SAR数据监测洪涝灾害的准确性。Geudtner等[3]、郭欣等[4]证明了利用InSAR技术提取洪水淹没区的可行性。目前，基于SAR影像的水体信息提取方法主要有阈值法[5]、结合DEM辅助提取[6]以及利用纹理信息提取[7]，其中阈值法简单易行，应用广泛。

面向对象提取方法是以具有相同特征的影像对象为基本单元进行处理的一种方法，通过对影像进行分割形成具有相同特征的对象，然后综合影像的光谱特征、几何特征以及各类指数等特征构建合理的知识库完成目标地物的提取，解决了传统提取方法中的“椒盐”噪声现象[8]。Laso等[9]、杨朝辉等[10]分别采用面向对象方法实现地物分类，证明了面向对象方法在精度上的优势。

本文采用Sentinel-1B SAR数据利用面向对象的阈值分割方法对2018年8月山东省寿光洪灾进行动态监测，并结合利用Sentienl-2A光学数据进行面向对象CART决策树分类所得到的研究区土地利用现状进行灾害评估，为灾后重建提供数据支持和参考。

1 研究数据和方法

2018年8月18-20日，受热带风暴“温比亚”的影响，山东省弥河流域遭受了特大暴雨，弥河上游的冶源水库、黑虎山水库和嵩山水库水位高涨，入库流量远大于出库流量，为确保水库安全，三座水库向下游泄洪，弥河及其支流所流经的县市受暴雨和上游水库泄洪的双重影响发生了严重的洪涝灾害，损失严重。本文选择弥河流域(36°04′～37°19′N，118°10′～119°10′E)作为研究区进行洪水监测，该研究区包括寿光市、青州市、临朐县和昌乐县(图1)。

图1 研究区地理位置及土地利用类型分类图Fig.1 Geographic location and land use types classification map of the study area

1.1 数据源

本文选取灾前、灾中、灾后3个时期的Sentinel-1B SAR数据进行洪水淹没面积监测，数据下载于欧洲航天局官方网站(http://scihub.copernicus.eu)。Sentinel-1B SAR数据采用Level-1级的地距多视(GRD)产品影像，该产品是经过多视处理并采用WGS84椭球投影至地距的聚焦数据[11]。数据获取时间和具体参数见表1。

表1 Sentinel-1B SAR影像Table 1 Sentinel-1B SAR image

本文选择平水期(2018年6月6日)的Sentinel-2A多光谱数据进行研究区土地利用类型分类，数据下载于欧洲航天局官方网站(http://scihub.copernicus.eu)。Sentinel-2A光学数据采用Level1-C级产品，该产品是经过几何精校正的正射影像。本文只使用了分辨率为10 m的2、3、4、8波段，数据具体参数见表2。

表2 Sentinel-2A光学影像Table 2 Sentinel-2A optical image

另外，本文还利用了研究区行政区划、DEM等数据来辅助进行淹没区的提取与分析。其中，DEM数据采用美国NASA和NIMA联合研制的30m分辨率的SRTM-1高程数据。

1.2 研究方法

1.2.1 数据预处理

本文使用欧空局发布的SNAP软件对Sentinel-1B SAR数据进行预处理，主要包括轨道纠正、辐射定标、几何校正、滤波、镶嵌、裁剪、后向散射系数转换等。几何校正采用SNAP内嵌的Range-Doppler算法，在校正的同时将影像重投影到UTM Zone50N/WGS84坐标系下，滤波采用Refined Lee算法，在抑制斑点噪声的同时保留地物的边缘信息。

Sentinel-2A的L1C级多光谱数据是经过几何精校正的正射影像，还需对其进行大气校正处理来得到L2A级数据。本文基于欧空局发布的专门用于生产L2A级数据的SNAP插件Sen2cor对Sentinel-2A数据进行大气校正。利用SNAP软件对大气校正后的影像进行地形校正，并同时重投影到UTM Zone50N/WGS84坐标系下。对Sentinel-2A的两景影像进行镶嵌、裁剪，得到2018年6月6日的高清卫星影像。

1.2.2 基于Sentinel-1B SAR数据的水体提取

在影像预处理的基础上，本文采用面向对象方法进行水体提取，主要包括影像分割和水体提取两个过程。本文使用eCognition Developer软件对研究区域进行面向对象水体提取，技术路线如图2所示。

图2 技术路线图Fig.2 Technology road map

本文采用eCognition Developer软件中的多尺度分割方法对影像进行分割。在进行多尺度分割时需要设置合适的分割参数：分割尺度、波段权重、光谱因子、形状因子、光滑度和紧致度。结合研究区水体分布情况，经过多次试验，最终确定出本文的分割参数：分割尺度为30，颜色因子为0.9，形状因子为0.1，紧致度为0.5，平滑度为0.5。在该参数下得到的分割效果最好，水体与非水体分离界限明显，图斑边界较完整。

水体在微波范围内较低的后向散射特性为在SAR影像上确定洪水范围提供了依据[12]，如：平静的水面在SAR影像上呈现黑色，而建筑物、植被等在SAR影像上则呈现明亮色调。为研究水体与其他地物在微波波段后向散射特性的区别，本文借助Google Earth高分辨率影像，在洪水未发生时(2018年8月8日)的SAR影像上均匀选择林地、耕地、住宅用地、设施农用地、水域5种占研究区比例较大的地物类型作为感兴趣区域，并统计了5种地物的后向散射特性，结果如图3所示。通过分析图3可知，水体的后向散射系数较其他地物的小。根据水体的这一特性，本文选择用影像的后向散射均值来构建知识库，然后利用阈值分割法将水体提取出来。

图3 研究区主要地物后向散射特性Fig.3 Backscattering coefficient of main lands in the study area

1.2.3 基于DEM数据的山体阴影剔除

由于SAR影像是侧视斜距成像，所以当雷达经过山区时会出现透视收缩、叠掩和阴影等现象。山体阴影的后向散射系数与水体相近[7]，在SAR影像也呈现黑色或暗色，在利用面向对象阈值法提取水体时无法完全将水体与山体阴影区分开来，所以提取的水体中混杂着大量的山体阴影，影响了水体提取结果的精度。

临朐县南边和青州市南边地形以山地为主，地势高低起伏，所以本文借助DEM数据通过对临朐县和青州市进行地形分析来剔除山体阴影的影响，主要思路为：使用ArcGIS对DEM进行处理，提取出坡度信息和高程信息，基于坡度和高程信息分别设置合适的阈值提取出山体阴影。

1.2.4 基于Sentinel-2A数据的土地利用类型分类

为分析研究区内各类地物的淹没情况，本文利用洪水未发生时(2018年6月6日)的Sentinel-2A光学影像对研究区进行土地利用类型分类。参考2017年国土资源部组织修订的国家标准《土地利用现状分类》(GB/T 21010-2017),并结合本文的研究目的和研究区实际情况，将研究区分为8类：耕地、林地、工矿仓储用地(工业用地、盐田)、住宅用地(城镇住宅用地、农村宅基地)、交通用地、水域(河流、湖泊、水库)、设施农用地(蔬菜大棚、花卉大棚)和裸土地。

本文利用eCognition Developer软件采用面向对象CART决策树方法进行研究区土地利用类型分类。首先对预处理后的Sentinel-2A数据进行多尺度分割，在多次试验的基础上，最终确定分割参数为：分割尺度为45，颜色因子为0.6，形状因子为0.4，紧致度为0.6，平滑度为0.4。在该分割参数下得到的分割影像区别明显且图斑完整，满足“最大同质性，最小异质性”的准则[13]，满足后续的分类要求。对照Google Earth高分辨率影像在eCognition Developer软件中选取各类地物的训练样本，然后选取Sentinel-2A影像的4个波段、归一化植被指数(NDVI)、归一化差异水体指数(NDWI)、长宽比(Length/Width)、形状指数(Shape index)、亮度值(Brightness)和同质性(GLCM Homogeneity)10个特征变量作为输入数据集，利用CART分类器对样本进行训练生成CART决策树分类规则，最后根据生成的分类规则对影像进行分类[14]。分类完成后再根据研究区实际情况通过人机交互解译的方法对分类结果进行修正,得到最终分类结果。

2 结果与分析

2.1 Sentinel-1B SAR水体提取结果

2.1.1 提取结果

通过分析影像的后向散射均值，确定出洪灾前(2018年8月8日)、洪灾中(2018年8月20日)、洪灾后(2018年9月1日)水体提取的阈值分别为-15.68 dB、-11.98 dB、-14.70 dB， 小于该阈值的为水体， 大于该阈值的为非水体，利用eCognition Developer软件中的Assign Class函数设置阈值范围将水体与非水体区分开来，然后利用软件中的平滑算法对提取结果进行处理。提取的水体结果中混杂着部分山体阴影，为减少山体阴影对水体提取结果的影响，本文结合SRTM-1 DEM数据对研究区进行高程和坡度分析，确定出地形起伏较大区域的高程阈值为302 m，坡度阈值为12°，通过阈值分割的方法提取出山体，进而剔除误提的山体阴影。将剔除山体阴影后的水体提取结果与原始影像叠加，如图4所示。

图4 寿光洪灾前、中、后水体提取结果轮廓示意图Fig.4 Outline drawing of water extraction results in pre-flood, flooding and after-flood of Shouguang

从图4可以看出，本文提取的水体信息与原始影像上的水体高度吻合，3大水库轮廓较完整，弥河流域除部分地方出现断流现象外，整体边缘信息较好。

2.1.2 精度评价

本文参照Google Earth高分辨率影像，通过人工目视解译的方法利用ENVI软件在研究区内随机选取了一定数量的样本作为精度评价的验证样本。然后利用eCognition Developer软件中的基于对象样本的混淆矩阵精度评价方法对水体提取结果进行精度评价，结果如表3所示。

表3 水体提取结果混淆矩阵Table 3 Confusion matrix of water extraction results

由表3可知，基于Sentinel-1B SAR数据的面向对象水体提取结果总体精度均在95%以上，Kappa系数均在0.85以上，水体提取结果比较理想，可以用于洪水灾害监测，这同时也表明基于SAR数据进行洪灾监测的有效性和可行性。提取结果中出现的误提和漏提现象主要有3方面的原因：1是受SAR成像方式的影响，研究区中的一些细小的河流、湖泊无法被提取出来；2是参考影像与Sentinel-1B影像获取日期间的微小差异造成的；3是DEM精度较低，利用DEM得到的山体阴影中混杂着水体，导致少量水体被当作山体阴影剔除。

2.2 Sentinel-2A土地利用类型分类结果

2.2.1 分类结果

研究区的土地利用类型分类结果如```图1所示，各类地物所占面积及比例如表4所示。分析分类结果可知，研究区内主要地物类型为耕地、林地、设施农用地以及住宅用地，分别占研究区总面积的37.3%，28.5%，11.2%和9.8%。

表4 研究区各类地物所占面积及比例Table 4 The area and proportion of different lands in the study area

2.2.2 精度评价

为验证面向对象CART决策树分类方法的准确性，本文基于Google Earth高分辨率影像在研究区内均匀选择了各类地物的验证样本，然后利用eCognition Developer软件中的基于对象样本的混淆矩阵精度评价方法对分类结果进行精度评价，结果如表5所示，包括总体分类精度、Kappa系数、混淆矩阵、用户精度和生产者精度。

表5 研究区土地利用类型分类结果混淆矩阵Table 5 Confusion matrix of classification results of land use types in the study area

分析表5可知，分类总体精度为85.7%，Kappa系数为0.835，分类结果与研究区实际情况基本吻合，可以用于后续的灾情分析。

2.3 洪水灾害评估

本文对研究区内各分县在洪水过程中的水体面积进行统计，结果如表6所示。监测结果显示，灾前水体总面积为305.98 km2，灾中水体总面积为504.57 km2，灾后水体总面积为371.08 km2。对比灾前灾中水体面积，寿光市水体面积增加最多，为167.44 km2，受灾状况最严重；临朐县、昌乐县、青州市水体增加面积依次为12.21 km2、10.29 km2、8.65 km2，这3个分县的水体面积增加相对较少，受灾状况较轻。对比灾中灾后水体面积，洪水发生12天后，昌乐县、临朐县和青州市的洪水基本消退，受灾最严重的寿光市地势低于其他3个分县，洪水短时间内难以完全消退。

表6 寿光洪水前、中、后研究区各分县水体面积Table 6 Water area of each county in pre-flood, flooding and after-flood of Shouguang

为研究洪水在研究区内的空间分布，本文利用ArcGIS软件对洪水过程中不同时间段的水体信息进行叠加分析，结果如图5所示。

分析图5(a)可知，研究区内水体面积增加的形式主要有河道拓宽、低洼地区积水和湖泊面积增加3种，受洪灾影响最严重的区域位于弥河两岸以及寿光市北部。结合图5(a)和图5(b)可知，临朐县、青州市和昌乐县三地洪水基本消退，水体面积下降较快。寿光市由于地势较低且弥河上游水库泄洪量巨大，洪水下降速度缓慢，退水周期较长。由图5(c)可知，寿光市是此次洪灾中受灾最严重的县市,尤其是北部沿海的原盐生产基地，受降雨和潮水水位的双重影响，出现严重洪涝灾害。另外，强降雨导致弥河两岸大量的蔬菜大棚被淹没，并且大棚内地下渗水现象严重，洪水消退缓慢。

图5 水体动态变化图Fig.5 Dynamic change map of water

本文结合洪水空间分布图与研究区土地利用类型分类图对研究区内各类地物的淹没情况进行分析，结果如表7所示。由表7可知，洪水淹没的主要地物类型为工矿仓储用地、交通用地和耕地，分别占总受灾面积的48.90%、13.90%、12.35%。被淹的工矿仓储用地主要集中在寿光市北部的大型盐场，淹没的交通用地集中分布在寿光市以及昌乐县东部地形平坦的地区，大多数为蔬菜大棚间的道路，研究区内淹没的耕地主要分布在寿光南部、青州东部以及昌乐北部，位于寿光市北部的大多数村庄受地势影响，许多房屋和养殖基地被淹没。

表7 研究区各类地物受淹面积Table 7 Submerged area of different lands in the study area

另外，如表7所示，受强降雨、海水涨潮和地下水位的多重影响，洪水过后，工矿仓储用地和交通用地仍有大量积水，淹没面积分别达38.06 km2、8.47 km2，结合洪水动态变化图可知，退水缓慢的工矿仓储用地主要为北部盐田区域，交通用地主要分布于被淹没的蔬菜大棚间。

分析不同时期的水体空间分布图可知，此次洪水主要分布在弥河沿岸以及弥河下游的寿光市。由于管理不善，弥河河道被大面积占用，房屋、耕地、养殖棚等违建随处可见，给弥河的正常泄洪造成隐患，所以政府应加大整治力度，拆除一切违法建筑，恢复河道原有宽度，提高河道泄洪能力。位于弥河下游的寿光市受到暴雨和上游泄洪的双重影响损失最严重，寿光市地形为半山区半平原，地势西南高东北低，所以在洪水到来之前位于地势低洼处的村庄应该及时撤离到附近地势较高的地区，尤其是弥河经过的村庄。另外，在雨季来临之前，各村庄应做好防洪工作，如修建大坝、蓄水池、引水沟渠等。在洪灾后应做好防疫消毒工作，采用无人机和喷洒车等方式持续对被淹没的村庄和农田等进行全方位消毒，避免灾后疫情的出现。

3 结论与讨论

根据水体的后向散射特性，本文基于Sentinel-1B SAR数据对2018年寿光洪灾进行遥感监测，统计出了研究区内各分县在洪水过程中的水体面积变化以及各地物的淹没情况，主要得到以下几点结论：

(1)SAR数据具有对地观测的全天时、全天候等特点，本文基于Sentinel-1B SAR数据得到的水体提取结果精度较高，总体精度均在95%以上，Kappa系数均在0.85以上，证明Sentine-1B SAR数据可以对受灾地区进行快速且准确的监测，在洪灾监测领域具有很大优势。

(2)研究区内水域面积总共增加198.59 km2，受灾最严重的区域位于弥河、水库、湖泊周围地势相对平坦低洼的地区，其中位于弥河下游的寿光市受淹面积最大，达167.44 km2。

(3)通过对研究区各地物类型受淹面积进行统计，工矿仓储用地、交通用地和耕地受灾情况最为严重，受淹面积分别为90.25 km2、25.65 km2和22.79 km2，其中被淹没的工矿仓储用地主要位于寿光市北部盐区。

Sentinel-1B SAR数据在洪水监测方面具有很大的优势，但单一的遥感数据无法保证监测结果的精度，所以在后期的研究中应该加强以下方面的研究：

(1)结合Sentinel-1A SAR数据，选择更多时相的数据对洪水过程进行监测，并充分利用Sentinel-1 SAR的双极化数据来增强水体特性，探讨更为快速准确的水体提取方法。

(2)在保证精度的前提下，探讨薄云覆盖下的光学数据水体提取方法，与SAR数据的监测结果相结合，在保证监测时效性的同时降低监测成本。