姜晗兵 邓文彬
(新疆大学,乌鲁木齐 830017)
河南郑州“7·20”特大暴雨洪涝灾害造成城市内涝、山体滑坡、河流漫堤等险情,导致重大人员伤亡及经济财产损失。2021 年7 月20 日8 时至21 日6 时,河南中北部地区出现历史极值降雨,郑州、新乡、开封、周口和焦作等地降雨强度较大。7 月21 日,河南省部分地区极端降雨形势持续,贾鲁河上游洪水严重超出其承载能力,部分地方溃堤。据评估,从郑州等地泄洪的洪水流量达到1 000~1 200 m3/s,而贾鲁河一般可承受的洪水流量仅为200 m3/s,周口市河段存在溃堤风险,态势非常危急。贾鲁河泄洪至周口市首先经过位于其西北部的扶沟县,因此本次特大暴雨及贾鲁河泄洪对扶沟县造成了一定影响。
遥感技术因其具有宏观、综合、快速、动态等特征,可为洪涝灾害监测预报提供非常重要的帮助,遥感技术在洪水范围提取方面已有非常多的应用[1]。洪涝灾害监测的关键在于水体范围的提取,遥感技术用于水体信息的提取主要分为光学遥感和雷达遥感,由于洪涝灾害发生的同时往往伴随着暴雨的出现,会有较厚的云层覆盖,因此难以获取高质量的光学影像[2-5]。雷达数据具有全天时、全天候、穿透性的特点,越来越多的学者将雷达数据运用在洪涝灾害的研究中[6-8]。阈值法、谱间分析法和多波段运算法是遥感上较常应用于水体范围提取的几种算法。Shu等[9]利用雷达影像得到洪水范围监测图,得出SAR 数据在洪涝灾害监测方面具有重要作用。李景刚等[10]利用欧洲环境卫星提供的数据提取了2007 年洞庭湖枯水期和洪水期两景影像上的水体范围,结果显示其改进的阈值法最优阈值的水体范围提取的精度高于双峰法和最大类间方差法。汤玲英等[4]基于Sentinel-1 影像,利用面向对象的方法绘制出西桂林会仙岩溶湿地的受灾前后洪水面积变化图。孙书腾等[11]将Sentinel-1 和Sentinel-2 数据相结合利用水体指数SDWI 提取出河南浚县“7·20”洪涝灾害前后的洪水面积,并对土地利用类别的灾情进行分析。
本文选取Sentinel-1 和Sentinel-2 影像数据,利用雷达和光学数据优势互补的特点,对研究区灾中、灾后洪涝灾害受灾范围进行快速提取,分析扶沟县在“7·20”特大暴雨影响下的受灾情况。
扶沟县隶属于周口市,坐落于周口市西北地区,总面积1 173 km2。截至2022 年,扶沟县共辖8 个镇、6 个乡;境内河流均属淮河流域,总长度206.4 km,年平均排涝量0.88 亿m3,最长河流为贾鲁河,总长为47.2 km,贾鲁河主要支流有双洎河、康沟河、卢义沟等。境内另有大浪沟、清水河、庙陵岗泓、东西黄水沟等;人工河流有幸福河、丰收河等。
在遭受极端强降雨及贾鲁河上游泄洪后,被誉为“蔬菜之乡”的扶沟县损失惨重,过水行政村达39 个,农田受灾面积达5.64 万hm2,受灾人口达56.1 万人,水利、交通、市政等基础设施受到不同程度影响,直接经济损失达28.4 亿元。
Sentinel-1是欧洲航天局哥白尼计划(GMES)中的对地观测卫星,由两颗卫星星座组成,有4 种成像模式,分别为SM、IW、EW、WV 模式,IW 模式是Sentinel-1 提供的主要数据[12]。Sentinel-2由两个星座组成,A星座于2015年6月开始使 用,B 星 座 于2017 年3 月 开 始 使 用。Sentinel-2 包 含13 个波段,分为10 m、20 m 和60 m 3 个空间分辨率等级的多光谱影像[13]。本文选用Sentinel-2数据是已经过正射校正和几何精校正的L1C级数据,选取其中3个空间分辨率都为10 m的B2、B3、B4波段,研究所需的Sentinel-1和Sentinel-2影像数据如表1所示。除此之外,还获取了扶沟县的行政边界数据及30 m分辨率的数字高程模型(DEM)。
表1 Sentinel-1和Sentinel-2影像数据
本文选取的是IW 拍摄模式下的level-1 级地距影像(GRD)数据,地距影像产品经过多视处理[14]。经判断,下载的影像数据都为升轨数据。本研究基于欧洲航天局提供的专业软件SNAP 8.0 对Sentinel-1 数据进行预处理,首先进行轨道校正,轨道校正时会自动下载精密的轨道数据并更新元文件中的卫星轨道状态信息,能够获取更准确的卫星轨道位置;再是辐射定标,这一步是为了将没有单位的后向散射信号转化为有单位的物理量,由于SAR数据具有云层的穿透性,因此辐射校正中只需进行辐射定标;而后进行相干斑滤波,在本次研究中运用的是Refined Lee 滤波器;之后进行影像裁剪,再地形校正,在地形校正这一步中不仅会给影像赋予真实的坐标信息,还会进行地形辐射校正,用的是30 m 分辨率的DEM 数据;最后进行分贝化,经过上述处理的后向散射系数的值较小,将雷达数据经过分贝化后,对雷达数据的分析会更加方便可用。
研究结果表明,在一系列对Sentinel-2 L1C 级产品进行大气校正的模型中,Sen2Cor 模型具有较高的精度。因此,本研究采用Sen2Sor 模型对三期Sentinel-2 L1C 产品数据进行预处理,利用SNAP 8.0 的重采样工具对B2、B3、B4这3 个波段进行重采样,输出为ENVI 格式,利用ENVI5.3进行波段合成、镶嵌、裁剪等操作。
本研究的主要流程分为3 个部分:研究区灾前土地利用分类、洪涝灾害淹没范围提取、洪涝灾害的空间分析。
利用支持向量机的方法对预处理后的Sentinel-2 影像进行土地利用分类。支持向量机分类(SVM)是一种建立在统计学习理论(SLT)基础上的机器学习方法[15]。支持向量机是一种适合小样本的学习方法,一般不会涉及概率测度及大数定律等,因此将分类和回归等问题变得简便了。支持向量的个数而不是样本的空间维数决定计算的复杂程度,这也在某种程度上避免了“维数灾难”的问题。支持向量机算法的学习过程可以概括为凸优化问题,因此能够运用已有的算法发现目标函数的全局最小值。参考《土地利用现状分类》[16],再结合扶沟县的真实土地利用状况,将其分为耕地、水体、建设用地和其他用地4 种土地利用类型。在4 月17 日的影像上随机选择一定数量的样本,判断每个类别样本之间的可分离性是否符合要求,可分离性均大于0.8,满足分类要求。利用ENVI 混淆矩阵工具评价分类精度,利用ROI 工具建立验证样本集,结果显示:总体的分类精度是95.85%,Kappa 系数是0.94,表明得到的灾前土地利用分类图可为后续的灾情分析 提供必要支持。2021 年4 月7 日Sentinel-2 土地利 用 分类图如图1 所示。
图1 2021年4月7日Sentinel-2土地利用分类图
由于传感器接收的微波后向散射的强弱不同,非水体表面一般比较粗糙,水体表面比较平整,因此在雷达影像上就会呈现出深浅不一的色调。水体表面较平整,因此微波的后向散射能力较小,呈现出黑色或者较暗的颜色[17]。
由于洪涝灾害一般具有突发性、短时性,因此对于洪涝灾害范围的提取应考虑到简单、快速。贾诗超等[18]根据雷达数据上水体和非水体的不同点及NDWI、NDVI两种指数计算方法的特点,探索出雷达影像的VV 极化数据结合VH 极化数据提取水体信息的初始关系式,通过持续的改进和实验逐渐完善提取水体信息的拟合关系式,最终探索出了基于Sentinel-1双极化数据的SDWI水体提取指数,公式如下:
式中:KSDWI为水体提取指数;VV、VH为Sentinel-1 极化影像。理论上讲,当KSDWI大于0 时为水体,KSDWI小于0 时为非水体,可根据实际情况,选择最合适的阈值。
式(1)的核心思想是通过将VV和VH极化数据相乘之后再乘以10,来扩大水陆之间的差异,然后再取自然对数,当自然对数大于1时,曲线的斜率也会慢慢变小,以此来找到合适的经验阈值。
经过预处理后的VV、VH极化影像在ENVI5.3中显示,利用波段计算器计算SDWI,从VV、VH、SDWI 3 个直方图的分布情况来看,SDWI 的双峰形态最为明显,VV 次之、VH 最差。因此,SDWI 最适合利用双峰法进行快速提取洪涝灾害的范围。
通过反复实验对比分析,最终确定提取灾中、灾后水体信息的阈值为-0.2。由于扶沟县为平原地区,没有山体,地形起伏不大,因此不考虑地形产生的阴影对水体信息提取的影响。通过决策树分类得到水体初步提取结果。
初步提取的水体范围会缺少空间连续性,需要进行小斑块的去除,首先进行聚类处理,聚类处理可以使原本的分类结果更加平滑,原理是利用数学形态学算子(腐蚀和膨胀),将邻近的类似分类区域聚类并进行合并。再进行最大值、最小值(Majority、Minority)分析,将较大类别中的虚假像元归到该类中,分析完成后一些较小的斑点被分到了背景类别中,使得水体范围提取结果更加平滑。经过以上处理,大量斑点像元被剔除,水体边缘得到了增强。灾中、灾后水体范围提取结果如图2和图3所示。
图2 7月27日水体范围提取结果
图3 8月8日水体范围提取结果
利用基于地面真实样本的混淆矩阵的方法对提取的灾中、灾后的水体范围进行精度评价,在Sentinel-2影像上随机生成200 个点作为精度验证的真实样本。随机点分布位置如图4 所示。由此基于SDWI 进行水体范围提取结果的精度为:7 月27 日的总体精度是97.6%,8 月8 日的总体精度是95.4%;灾中、灾后提取的水体范围精度均高于95%,因此满足实验需求。
图4 随机点分布位置图
将洪涝灾害前基于Sentinel-2 影像提取出的土地利用图及灾中、灾后提取出的水体面积统一投影在WGS_1984_UTM_Zone_50N坐标系下,得出研究区灾前、灾中、灾后的水体面积变化情况,结果如表2所示。
表2 研究区灾前、灾中、灾后水体面积变化情况km2
结果表明,研究区在2021 年4 月17 日、7 月27 日、8 月8 日的水体总面积分别为4.068 km2、36.468 km2、18.770 km2,可以看出灾中水体面积比灾前的水体面积扩大了近9 倍。从7 月20 日暴雨开始,叠加贾鲁河上游的泄洪,到7月27日水体面积急速增长到36.468 km2,洪涝灾害面积高达32.382 km2,随着强降雨的过境和泄洪结束,洪水的面积也在逐渐减少,8 月8 日水体面积缩减到18.770 km2。对比灾前、灾中、灾后各乡(镇)的水体面积,曹里乡受灾面积最大,水体最大变化面积达到12.630 km2。扶沟县地处贾鲁河与双洎河交汇处,受强降雨和两条河流上游泄洪影响,扶沟县发生严重洪涝灾害,其中曹里乡受灾最重,26个行政村全部受灾,11个行政村受灾较为严重,房屋进水,受灾人口42 667人,直接经济损失5.95亿元;白潭镇受灾情况次之,水体最大变化面积6.432 km2,由灾中、灾后水体面积变化可知,白潭镇洪涝灾害面积消退较快,8月8日基本同灾前面积接近;练寺镇、汴岗镇、韭园镇、城郊乡的受灾情况较轻,最大水体面积变化在2~4 km2,但这几个乡(镇)灾情较为持久,洪水消退较慢;其余几个乡(镇)由于距离贾鲁河和双洎河较远,因此受此次暴雨和上游泄洪的影响不大。
扶沟县本次洪涝灾害事发突然,7 月19—22 日,持续超强降雨造成的汛情、洪峰均达到扶沟县气象、水文有记录以来的最高值,周口市水利局2021 年7 月21 日2 时起,将水旱灾害防御Ⅳ级应急响应提升至Ⅱ级。贾鲁河上游向下游加大泄洪量,同时南水北调干渠水量也向贾鲁河下游泄洪,周口市贾鲁河出现近40 年来的最大洪水。曹里乡位于贾鲁河与双洎河交界处,受此次泄洪影响最大,从Sentinel-2 影像上可以看出,贾鲁河及双洎河在曹里乡境内河堤发生了不同程度的溃堤,导致曹里乡的行政村全部受灾。基于灾前土地利用分类及灾中水体范围提取结果,统计出全县农田受灾面积为22.793 km2,主要集中在曹里乡。
本文选择周口扶沟县为研究区,基于Sentinel-1 和Sentinel-2 雷达影像快速提取研究区受洪涝灾害影响范围,本研究对洪涝灾害发生时快速获取受灾情况有一定的借鉴意义。主要结论如下。
(1)基于雷达数据进行洪涝灾害范围提取的研究中,阈值分割法提取速度较快,特别对于平原地区,无山地阴影影响的区域,可以快速、准确地提取出洪涝灾害的范围,在洪涝灾害应急监测方面发挥着非常重要的作用。
(2)本文利用Sentinel-2 光学影像结合支持向量机算法进行灾前的土地利用分类,精度达95.85%。利用Sentinel-1雷达影像结合SDWI对灾中、灾后的水体范围进行提取。最终提取结果的总体精度分别为97.6%和95.4%,均在95%以上。表明利用Sentinel-1雷达影像结合SDWI能够快速准确进行洪涝灾害范围的提取。
(3)扶沟县灾中的水体面积呈现出“突增”的变化趋势,灾中水体面积相比于灾前增加了32.382 km2,灾中水体面积为灾前水体面积的近9倍。
(4)扶沟县境内两条河流贾鲁河、双洎河在曹里乡汇流,由于连日强降雨及贾鲁河上游泄洪,部分河堤河水漫溢决堤,导致曹里乡受灾最严重,水体最大变化面积达12.630 km2。
运用主被动遥感相结合的方法在洪涝灾害发生时进行水体范围的提取,这在灾害发生前后对灾情进行评估及制定救援策略方面发挥着非常重要的作用。洪涝灾害发生时一般伴随着特大暴雨,由于此时云层较厚、水汽较为充足,仅利用光学影像进行水体范围提取会受到一定的限制。而具有穿透性、不受天气影响的雷达影像此时能发挥更好的作用。但由于雷达影像的成像原理比较复杂,数据处理的难度较大,因此使用的领域和范围受到限制,所以,本文采用雷达和光学遥感协同处理的方法进行研究区洪涝灾害范围的提取分析。
影响本研究精度的主要因素有:利用阈值分割法进行水体信息提取时,阈值比较难确定,需要进行反复试验:利用光学影像对雷达影像提取的水体范围进行精度验证时,验证样本与提取影像存在时间差。因此,未来可以加强以下几方面的研究。
(1)增加灾害发生过程中可用的遥感数据源,由此可以解决单一遥感数据在洪涝灾害监测方面的不足,缩短洪涝灾害监测的时间间隔,更好地绘制洪涝灾害期间淹没情况。
(2)加强不同场景水体散射特性的研究,对有些受灾情况不是很严重的农田范围进行提取,精确洪涝灾害淹没范围。
(3)加强多源遥感数据融合的研究,利用光学遥感和雷达遥感优势互补的特点提高洪涝灾害范围提取的准确性。