李勇永,郭 涛,张晓萍,王学恭
(1. 洛阳师范学院 国土与旅游学院,河南 洛阳 471934;2. 四川省农业科学院 遥感与数字农业研究所,四川 成都 610066;3. 天水师范学院 资源与环境工程学院,甘肃 天水 741001)
Google Earth Engine(GEE)是一个全球尺度地理空间数据云服务平台,它继承了多种地理空间数据并提供了超过1 000种的数据类型和算子,用户无需下载数据即可直接在线分析。近年来,众多学者借助遥感云计算平台GEE 来获取水体[1-10]和生态环境信息[11-13]。目前,基于GEE平台评价洪灾对郑州生态影响的研究较少,洪灾对人类造成了巨大伤害,对不同环境下植被生长的影响如何?在GEE平台高效的洪灾动态变化监测基础上开展灾后生态影响评估显得尤为重要。
本文综合多源、多时相遥感数据,从不同角度提取郑州市7·20洪灾水体信息,并对洪灾影响下的城市生态环境变化进行监测与评估,进而表明多光谱和雷达影像数据都能在洪灾监测和生态破坏评价中发挥重要作用,为决策者提供支持,让政府或公众可以直观了解洪灾的严重程度及其对生态的影响程度。
郑州市地处112°42′~114°14′E,34°16′~34°58′N,位于河南省中部偏北地区,地势西南高东北低(图1)。全市总面积约7 446 km2。郑州境内以黄河、淮河水系为主,旱涝灾害频繁,气候变异性明显。2021年7 月20 日—7 月30 日,河南省遭遇长时间大规模降雨,郑州“7·20”特大暴雨强度和范围突破历史记录,全市城乡大面积受淹,城镇街道洼地积涝严重、河流水库洪水短时猛涨、山丘区溪流沟道大量壅水,形成重大自然灾害[14]。
图1 研究区域地形图
本文选用Senitnel-1、Sentinel-2 以及MODIS 影像数据提取郑州市洪水灾害区域并进行灾后生态评价。选用欧盟委员会联合研究中心(JRC)制作的全球地表水月尺度产品[15]与洪水灾发生前、中、后的洪水提取面积进行空间分布对比。选用武汉大学黄昕教授团队制作的近30 a中国土地覆盖数据产品(CLCD)[16]分析不同土地覆盖条件下洪水灾害对区域生态环境的影响。为便于研究,将2021 年6 月6 日—2021 年7 月20日、2021年7月20日—2021年7月30日、2021年7月30日—2021年9月30日界定为洪水发生前、中和后期3个阶段。文中涉及数据集如表1所示。
表1 遥感数据集
洪水淹没导致地表覆盖类型发生突变,改进的归一化差异水体指数(MNDWI)[17]可用于目标时段的水体指数计算,公式如下:
MNDWI=(GREEN-MIR)/(GREEN+MIR) (1)式中,GREEN为绿光波段;MIR为中红外波段;分别对应于哨兵2 号多光谱数据的波段3 和波段12。该方法对于植被、建筑物及裸露地表均具有较好的抑制作用[18],但是由于洪水降雨期间会出现云雾遮挡等情况,多光谱成像易受影响。
洪水发生时,SAR影像中地表面散射体的雷达后向散射系数在时序上会呈现异常波动,影像中亮度代表回波强度,水体回波强度较小,陆地回波强度较大,水体在影像中呈现暗色或黑色,陆地呈现灰白色或黑灰色[19]。SAR 影像虽然具有较强的穿透性,能克服云雾干扰,但受地形因素影响较大,往往会把地形阴影识别为水体。
水体与陆地其他地物之间的差异在MNDWI 指数和雷达后向散射值的灰度图像中较为突出,但还需使用阈值对图像进行分割。基于经预处理的影像灰度直方图,以目标和背景间方差最大为基础假设,采用最大类间方差法动态确定图像的全局最优阈值,实现水体提取[7]。
从绿度、湿度、热度及干度四个方面建立洪水灾区的遥感生态质量评估体系,并采用主成分分析法构建综合遥感生态指数[20-22],即:
式中,RESI 为遥感生态指数;NDVI 为植被指数;WET 为湿度;LST 为地表温度;NDBSI 为土壤指数。文中采用的主成分分析是一种采取依次垂直旋转坐标轴的方法,将多个变量的信息通过线性变换集中到少数几个特征分量的多维数据压缩技术[23]。该方法的最大优点是集成过程中,各分量指标的权重主要根据数据本身的性质以及各指标的贡献度客观确定,避免人为确定,造成结果偏差[24]。
植被覆盖率(FVC)指包括地上叶、茎、枝在内的植被垂直投影面积占统计区域总面积的百分比,是衡量植被生长状况的重要指标[25-26]。获取区域植被覆盖率及其变化,对揭示地表生态环境空间变化规律、探索变化驱动因素及评价区域生态环境具有重要意义。文中基于NDVI,通过像元二分模型来估算不同时期的植被覆盖率。
式中,NDVIsoil为纯土壤像元的植被指数;NDVIvege为纯植被像元的植被指数。
基于GEE平台对洪水灾害发生前、中、后的水体进行提取,并采用遥感生态指数和植被覆盖率评价洪水对区域生态环境的影响,关键过程如下。
1)在GEE 平台中筛选出研究区洪水灾害发生前(63景)、中(14景)、后(89景)的Sentinel-2 MSI数据,然后以无云像元重构合成3 个目标时段各自的最小云量影像集合,并求取集合中每幅合成影像的MNDWI 值,最后基于MNDWI 指数集合求取MNDWI均值影像,以克服云雾遮挡。
2)在GEE 平台中筛选出研究区洪水灾害发生前(13 景)、中(3 景)、后(20 景)的Sentinel-1 SAR GRD 数据,然后选取VV 极化波段并掩膜掉边缘无效值形成影像集合,最后基于经过掩膜的VV 波段集合求取SAR均值影像,以克服地形影响。
3)采用最大类间方差法确定MNDWI均值影像和SAR均值影像的水陆分割阈值,分别实现水体信息粗提取,然后联合SAR影像和多光谱影像处理结果,当两个提取结果在同一地理位置都是水体时,则判定为水体,否则都为非水体,以获取水体信息的精细提取结果。
4)采用吉林一号高分影像,基于人工标注真值进行定量验证;采用JRC水体数据集,基于空间分布进行定性验证;采用GIS 空间统计技术分析洪水淹没情况和不同土地覆盖类型的淹没程度。
5)在GEE 平台中筛选出研究区2019—2022 每年8 月份的MODIS 数据产品,分别用于计算热度(MOD11A2 每年各4 景),绿度(MOD13A1 每年各2景)、 湿 度(MOD09A1 每 年 各4 景) 和 干 度(MOD09A1每年各4景),然后通过像元筛选合成最小云量影像以克服云雾遮挡,通过JRC水体数据进行掩膜以消除大量水体对遥感指数主成分计算的影响,最后基于研究区无云合成影像计算绿度、湿度、热度和干度,并用主成分分析计算获取研究区2019—2022年各年相同时段的遥感生态指数,将2019、2020 和2022三年的均值作为非受灾年遥感生态指数。
6)在GEE 平台中筛选出研究区2019—2022 每年8—9 月份的Sentinel-2 MSI 数据(分别有89、84、84、76 景),然后通过像元筛选合成每年的最小云量影像以克服云雾遮挡,并基于合成的最小云量影像计算NDVI,最后通过像元二分模型来估算2019—2022年各年相同时段的植被覆盖率,并将2019、2020 和2022三年的均值作为非受灾年植被覆盖率。
7)利用同期对比法评价灾后不同土地覆盖类型的生态状况。
3.1.1 洪水监测结果验证
利用吉林一号高分影像的人工标注真值进行精度定量评价,其中水体和非水体标注各120 个(图2)。通过建立混淆矩阵计算出评价洪水监测结果好坏的总体精度及表征水体提取一致性程度的Kappa 系数,对研究结果进行验证评价(表2)。
图2 吉林一号卫星影像与人工标注真值(Copyright,2021年7月24日由长光卫星技术有限公司提供)
由表2 可知,提取结果达到高度一致的水平,存在个别错分的原因有:①洪水期间地表覆盖的形态发生剧烈变化,影像数据集无法精确识别;②水体提取是以洪水发生期间,无云像元中值合成目标时段的最小云量影像为基础,所得结果表征的是整个洪水期间地表水体覆盖基本情况,而以拍摄于7月24日的高分影像进行验证会客观地导致部分验证误差。
表2 监测结果验证情况
将JRC的地表水空间分布与洪水发生前中后各阶段的地表水空间分布进行提取精度定性评价(图3)。经空间统计可知,JRC 的地表水面积为196.27 km2,洪水发生前为148.64 km2,洪水发生时为976.66 km2,洪水发生后为148.3 km2。JRC 计算的水域面积为全年总的水流量与季节性水的总和,文中计算的水域面积仅为6 月1 日—7 月20 日日的地表水面积,其总量小于JRC 的面积。通过自西向东统计不同经度范围内的面积(图3)可知,采用方法与JRC 保持了较高的空间一致性。
图3 洪灾前、中、后与JRC的水域面积对比
3.1.2 洪水淹没情况分析
为准确了解不同土地覆盖的受灾状况,从洪灾期间水域部分去除JRC水域,得到实际的洪水覆盖区域,并分析该部分水域内的土地覆盖状况(图4和表3)。结果表明,主要受灾的土地覆盖类型为耕地和不透水面,其中,耕地淹没面积为386 km2,占区域总面积的比例最大,不透水面淹没面积达274 km2,受影响的比例最大。
表3 主要土地覆盖受灾情况统计
图4 洪灾前、中、后与JRC的水域面积对比
3.2.1 遥感生态指数变化
基于空间化的绿度、湿度、热度和干度等遥感生态指数(图5),通过主成分分析获取2019—2022 年每年8 月的RESI,然后将非受灾年的RESI 均值与受灾年2021 年进行对比,以综合反映研究区在受灾后的生态状况,将RSEI 以0.2 为间隔分为I(0~0.2)、Ⅱ(0.2~0.4)、Ⅲ(0.4~0.6)、Ⅳ(0.6~0.8)、V(0.8~1.0),分别代表生态状况的差、较差、中等、良、优5个等级(图7)。由非受灾年至2021年研究区同期遥感生态指数的等级演变统计结果(表4)可知,2021年8月,生态状况等级提高的面积为2 060 km2,占比27.5%;等级退化的面积为422 km2,占比5.6%;等级不变的面积为4 616 km2,占比60.9%。从遥感生态指数的角度看,研究区生态状况总体较为稳定,呈转好趋势。以耕地为例,研究区耕地占比60%以上,大部分位于东部地区,但淹没比例较小,一方面降雨过后耕地湿度变大,另一方面转好的墒情促使以夏玉米为代表的作物长势向好,耕作区植被覆盖率变大。以上两方面对生态环境质量起正反馈作用,能够促进生态环境的好转[27]。
表4 非受灾年至2021年洪水灾区遥感生态指数等级转移矩阵/km2
图5 主要土地覆盖受灾分布图
图6 非受灾年与2021年的遥感生态指数各项指标变化
图7 洪水灾区遥感生态指数变化
3.2.2 植被覆盖率变化
对时空范围内经过预处理的哨兵2 号多光谱影像进行植被覆盖率(FVC)提取,其中包含时间筛选、空间筛选、去云和调用NDVI 函数,最后以各时段内的NDVI 均值合成结果为基础,评价地表植被覆盖情况。将FVC 以0.2为间隔分为I(0~0.2)、Ⅱ(0.2~0.4)、Ⅲ(0.4~0.46)、Ⅳ(0.6~0.8)、V(0.8~1.0),分别代表生态状况的差、较差、中等、良、优5个等级(图8)。由非受灾年至2021 年研究区同期植被覆盖率等级演变统计结果(表5)可知,生态状况等级提高的面积为1 704 km2,占研究区总面积的23.8%;等级退化的面积为1 725 km2,占研究区总面积的24.1%;等级不变的面积为3 730 km2,占研究区总面积的52.1%。从植被覆盖率角度看,灾后生态环境质量总体稳定。由良转为优(700 km2)在等级提高的类型中贡献较大,表明洪水过后植被覆盖率较高的区域生态状况进一步改善,这是因为洪水既能为植被提供水分,又能将养分带上来在一定程度上促进植被生长。
表5 非受灾年至2021年洪水灾区植被覆盖率等级转移矩阵/km2
图8 洪水灾区植被覆盖率变化
3.2.3 洪灾影响下不同土地覆盖的生态变化
不同土地覆盖对洪水灾害的敏感性差异较大[28],为获取各覆盖区域生态对洪水的响应,首先排除JRC水体区域,然后统计淹没区与非淹没区不同土地覆盖的遥感生态指数和植被覆盖率(表6、7)。
遥感生态指数区域统计(表6)显示,淹没区和非淹没区的生态总体状况都趋向好转,均值增长在11%以上。其中,耕地和不透水面的增长最为显著,林地和灌木覆盖条件下的生态变化不明显。
表6 不同土地覆盖的遥感生态指数变化/%
植被覆盖率区域统计(表7)显示,淹没区和非淹没区的生态状况总体趋向好转,淹没区均值增长24%,非淹没区均值增长12%。耕地和不透水面的生态状况显著转好;裸地区域没有植被覆盖;林地、灌木和草地的生态状况总体下降,尤其是草地类型,在淹没区下降60%,在非淹没区下降22%。需要说明的是,表7 中不透水面区域的植被覆盖率虽然在洪水过后都显著上升,但都在30%以下,根据不透水面与植被覆盖率在城市建成区呈负相关关系,植被覆盖率均值低于30%可划分为不透水面类型[29],证明研究结果具有较高的可信度。
表7 不同土地覆盖的植被覆盖率变化/%
总体来看,受人类活动影响较大的耕地和不透水面覆盖区的生态状况均有好转,受人类活动影响
较小的林地、灌木、草地覆盖区生态状况轻微变差;基于较高分辨率哨兵2 号多光谱影像获取的植被覆盖率对不同土地覆盖下生态变化的辨别能力要优于基于中分辨率MODIS 数据获取的遥感生态指数。
基于GEE 云平台,使用MNDWI 函数并参考最大类间方差法对Senitnel-1和Sentinel-2影像数据进行处理,提取洪水灾区的水体信息,然后采用MODIS 和Sentinel-2 MSI 数据分别计算非受灾年与2021 年灾后同期的遥感生态指数和植被覆盖率,实现生态状况快速监测分析。主要结论如下:
1)采 用GEE 联 合Sentinel-1 SAR GRD 和Senti⁃nel-2 MS数据集获取灾区的水体覆盖精度较高;耕地受灾面积最大,不透水面受灾比例最高。
2)相比非受灾年,研究区生态状况变化整体趋势稳定,灾后生态状况等级提高部分的占区域总面积的20%以上。
3)淹没区和非淹没区的生态状况都趋于好转,与人类活动密切关联的耕地、不透水面覆盖区域的生态状况呈显著好转,受人类活动影响较小的林地、灌木及草地的生态状况呈轻微恶化。