周 帆1,2张文君1 雷莉萍2 张智杰2 郭晓东 汪晓帆

(1.西南科技大学环境与资源学院 四川 绵阳 621010)(2.中国科学院遥感与数字地球研究所数字地球重点实验室 北京 100094)(3.中国土地勘测规划院国土资源部土地利用重点实验室 北京 100035)

引言

目前，洪水灾害被认为是世界上危害性最大的自然灾害之一，其造成的损失在全球所有自然灾害中占40%。洪水往往发生在人口密度高，河流和湖泊集中，降雨量大，农业种植的程度高的地方[1]。其对农业的影响主要表现为淹没耕地的大量的积水短时间内难以褪去，农田的受灾面积和灾后恢复状况均随着洪灾发生和褪去的时间不同而反应出不同的动态特征。这些特征是对灾情动态监测以及灾后损失评估的重要依据。

由于洪灾常伴随着持续的强降雨，具有突发性而不易实地调查。传统的调查和监测方法无法取得实时可靠的数据。遥感技术因其大范围、全天候全天时的特点被广泛应用于洪灾的监测，利用多时相的遥感影像可以实时监测洪灾的发生状况。喻光明等人分析了四湖流域“91.7”暴雨分区水量平衡及地面水情实况,利用遥感技术研究了该区“91.7”暴雨的洪水势态,包括洪水淹没范围估算及其时空特性分析,据此作出了灾情评价[2]。张国平等人利用遥感,对多种地表水体如河流和湖泊进行空间识别、定位及定量计算面积等研究内容开展洪水水情监测[3]。王云秀等人利用高分与资源卫星对暴雨后河北省中南部水库面积进行了监测分析[4]。李健等人利用利用FY-3A、HJ-1A/B和EOS多源卫星遥感数据,结合地面气象观测数据和基础地理信息数据,对2010年7月下旬至8月初吉林省特大暴雨导致的洪涝灾害进行了监测[5]。但在现实中洪灾对区域产生的影响不仅有农田淹没范围信息也应分析其受灾发生的时间、持续时间以及灾后恢复的动态信息。

本文利用多源遥感卫星的多时相数据对洪灾发生的时间、范围、不同区域的受灾程度以及灾后恢复状况进行动态监测。通过时序变化检测提取农田受灾动态信息的方法，为受灾状况分析以及灾后恢复评估提供决策的重要科学依据。

一、数据采集

由于2017年5月25日开始的季风性降雨导致了斯里兰卡境内山洪暴发，这次水灾是自2004年印度洋海啸以来，斯里兰卡遭遇的最大一次自然灾害。我们以洪水暴发的集中区域斯里兰卡南部马塔勒区作为研究区域(图1)。

图1 研究区域的地理位置和卫星遥感彩色影像

(2016年1月29日Landsat-8卫星观测的可见光和红外光(OLI6-5-3波段)波段合成的影像)

本文收集了2017年和未受灾的前一年2016年MODIS观测数据处理的16天合成NDVI数据产品(表1)。由于研究区域处于不同的轨道，将数据拼接后以16天为周期的范围内由受云影响最小以及距离星下最近的最佳像元合成，每年有23个时相。该数据产品来源于美国NASA的EOS系列卫星Terra观测获取的植被指数产品，全称为MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid(MOD13Q1)[6]。本文收集的NDVI数据产品是该植被指数产品之一，数据的空间分辨率为250m。16天合成的NDVI数据是由最佳像元选择合成算法处理得到的数据产品[7-8]。

作为实验的验证分析，还收集了灾中和灾后时期的GF-3和Sentinel-1观测数据以及灾期Sentinel-2光学观测数据和灾前Landsat影像数据。

GF-3数据来自于由中国高分应用技术中心分发的成像方式为精细条带1(Fine Strip I,FSI)模式，极化方式为双极化的Level-1A级数据。

Sentinel-1和Sentinel-2数据产品来自于欧空局(European Space Agency，ESA)官方网站：(https://sci-hub.Copernicus.eu/dhus//home)的Level-1地距数据(Ground Range Detected,GRD)，成像方式为干涉宽幅(Interferometric Wide swath,IW)模式，极化方式为VV的雷达数据。Sentinel-2的数据为经辐射校正和正射校正处理后的Level-1C级数据产品。

Landsat-8 OLI(Operational Land Imager)光学数据，来自于美国地质勘查局(U.S.Geological Survey-USGS)发布的经辐射校正和地形几何校正处理后的L1T数据产品。

表1 多颗卫星遥感观测数据产品信息

二、洪灾动态信息提取方法

本研究根据上述由研究区域样点的NDVI时序，本文提出了通过由洪灾导致的NDVI的变化特征来监测和提取耕地受灾信息的方法。其主要步骤包括对最初合成产品中异常点的剔除以及插值处理、多时相NDVI变化特征的动态监测、受灾程度以及灾后恢复状况的评价。

图2 基于多源卫星遥感洪灾动态信息提取方法流程

(一)NDVI时序变化特征

斯里兰卡南部主要栽培季节Maha，从9月到3月。主要灌溉的次级Yela季节从4月到9月，主要种植作物为水稻。以马塔勒地区2017年受灾区域为例，根据当地灾害管理中心(DMC)数据报告从中选出6个不同特征的农田样本区域并提取2017受灾年和2016年NDVI多时相数据。对比两年的NDVI变化情况(图3)。

图3 受灾2017年(红色)和未受灾2016年(蓝色)样点NDVI时序特征

灾前的持续强降水使土壤水分过度饱和以及灾后大面积农田被水淹没，严重影响农田作物的生长，该变化会通过遥感卫星监测数据中的地表反射率随时间序列的波动表现出来，具体表现为NDVI值随洪灾持续降低。因此根据农田的受灾变化与NDVI变化的关系特征可提取洪灾的动态信息。

1.对比未受灾年NDVI值，不同样区受灾年NDVI值开始降低的时期基本一致

根据灾害中心(DMC)的报告得知洪灾发生的日期为5月25日。如上图(a)(b)(c)(d)所示，受灾年NDVI值在5月末6月初下降趋势时分明显，表明洪灾发生后立刻对农田作物生长产生了严重的影响。

2.不同区域受灾后NDVI值降低的持续时间不同

将两年同时相的NDVI做差值对比发现受灾年同比减小约0.18，统计连续大于此阈值的时相数，其中(a)(b)(c)(d)分别呈现了8、7、4、3个时相，其中(c)(d)的NDVI值在受洪灾的影响降低，灾后又回升到了2016年的NDVI值水平。图1(a)(b)连续有9和8个时相的NDVI差值连续大于阈值，灾后一直持续降低至年末没有回升，该类区域表明在洪灾影响过后农民由于严重的破坏而放弃后续的耕作。由此可见，连续大于阈值的时相数可以反映受灾的严重程度以及灾后的恢复状况。

3.NDVI值异常降低与升高的区域

除上述受灾年NDVI值随洪灾降低的区域之外，也存在NDVI值异常降低或升高的情况。如图(e)从年初就一直低于为受灾年或一直高于为受灾年的变化，其因素可能为农田变为整年休耕或用于建筑用地。如图(f)所示异常升高又突然降低的区域可能为整个时相周期内均为多云天气导致。

(二)灾害动态信息提取

将上述耕地洪灾动态信息提取方法应用于斯里兰卡马塔勒地区受到洪灾2017年和未受灾的2016年MODIS多时相NDVI数据进行了方法验证。首先为保证数据的准确性对收集到的研究区域两年23个时相的MODIS数据进行异常点剔除以及线性插值等预处理；根据2017年与2016年NDVI的差值变化提取马塔勒地区耕地洪灾动态信息。

对比两年NDVI灾后同时相的变化差值数据表明，2017年受灾后的NDVI值比2016年同时相总体相差0.18。以此0.18为阈值，计算出NDVI降低时间Dt以及连续低于阈值的时相数Dn，根据这些特征可以评估作物受暴雨降水灾害的影响程度以及受灾后作物恢复情况。MODIS-16天NDVI产品数据虽然由最佳像元合成方法处理得到，很大程度上减小了云的影响[9]，但由于斯里兰卡的热带季风气候，雨季为每年5月至8月和11月至次年2月，持续时间很长。因此会出现在合成周期内均为多云天气的状况无法反应样点的真实数据，从而影响到NDVI样点的正常变化规律。因此在提取这些特征值之前，首先对MODIS-NDVI数据存在的部分异常进行剔除和插补。

在NDVI多时相变化中，按照农作物的生长到成熟收割，NDVI的变化呈现的是逐渐增大然后减小的变化规律[10]。NDVI的变化会随着作物在收获季或灾害影响而减小但不会在后一时相出现突然增大的现象。因此，若样点中连续时相之间出现突然增大和减小且幅度超过0.18的情况，则将其判定为受云层影响而产生的异常值并进行剔除。同时结合异常点前后时相NDVI值进行线性插值，由此减小异常值对多时相变化特征分析的影响。

将处理后得到2016年和2017年分别23个时相的NDVI数据集。通过总时相为N的受灾年(NDVIs[i])与未受灾年(NDVIc[i])相同区域与时相(i)中连续两个时相差值ΔNDVI[i](式(1)，i=1,2,…,N)大于阈值(T)作为受灾开始时间Dt。并将其受灾后持续时相数Dn(式(2)-式(4))作为受灾持续时长。

(1)

TNDVI[i]=

(2)

Dt=i,TNDVI[i]=1 and TNDVI[i+1]=1 i=1,2,…,N-1

(3)

Dn=Num(TNDVI[i]=1 Until(TNDVI[i]=0))i=Dt,Dt+1,…,N

(4)

根据上述采集到的样点区域NDVI随时间变化的特征结合通过阈值得出的受灾开始时间Dt以及受灾持续时间Dn，可将整个区域受灾情况划分为4种类型。

1.灾后恢复型：受灾开始时间Dt与洪灾发生时间基本一致。同时，受灾持续时间Dn小于60天的样点。

2.灾后弃耕型：受灾开始时间Dt与洪灾发生时间基本一致。同时，受灾持续时间Dn大于60天的样点。

3.非灾害型：研究区域中2017年NDVI值的减小除灾害原因外，还有其他变化的因素如耕地变为建筑用地、休耕地等类型。根据斯里兰卡的雨季变化规律，只有在5月(即i=9)才有可能发生持续强降水而引发洪灾的规律，以Dt必大于9或者Dn必小于(N-9)为条件限制，可以判定如图所示的NDVI值变化是由非灾害所引起的。

(三)结果与分析

分别为按上章所述方法逐像元处理，并按照像元分辨率得到不同特征像元的面积总数统计得到的受灾开始时间和持续天数的结果。图4的开始时间是根据NDVI数据的每个时相间隔16天，将开始时相Dt进行日期换算后的结果。图5的持续天数也是根据每个时相间隔16天由时相个数Dn换算的天数(图4)。

图4 研究区受灾信息直方图：(a)开始时期、(b)持续天数

同时，根据NDVI的变化特征以及受灾开始时间Dt、受灾持续时间Dn，根据上述不同受灾类型的划分标准，对整个受灾区域的进行分类提取结果(图5)。

图5 研究区域受灾类型提取统计图

1.从受灾开始时间提取结果可以看出，大部分受灾区域的开始时间在5月24日至6月9日之间，这也与灾害中心(DMC)报道的洪灾发生时间5月25日基本一致。部分区域洪灾发生前开始降低的现象可能是灾害来临之前的持续降水已经对部分区域的作物生长产生了影像。

2.从受灾持续时间提取结果可以看出，大部分受灾区域影像持续时间在49至65天之间。由此可见，洪灾对大部分区域作物的影响持续了两个月或更久，最久的区域持续了6个月。通过受灾持续时间的统计可以直观的反应不同区域受灾的严重程度，持续时间越长则表明受灾程度越严重。该数据可以为救灾计划实施提供有力依据。

3.从受灾类型的提取结果可以看出虽有48%的受灾区域在灾后恢复到了为受灾年的水平，但仍有38%的受灾区域由于影响较为严重在灾后很长一段时间内没有恢复到年前的状态，其中10%的非灾害类型根据联合国粮食及农业组织(Food and Agriculture Organization of the United Nations，FAO)出版的斯里兰卡2017作物和粮食评估报告中提到由于多年的干旱导致斯里兰卡粮食大幅减产使得部分农田被开发为建筑用地，在高海拔地区则通过休耕改种茶树等市场价值更高的作物。通过划分不同的受灾类型可以直观的反应不同地区灾后恢复情况[10]。该数据可以为灾后损失评估提供帮助。

三、对比与验证

SAR由于具有全天候全天时的特点，在极端天气条件下也能够在第一时间获取灾区的影像信息，为救灾工作以及灾情的预警评估提供重要的依据[11]。本研究利用灾中和灾后GF-3和Sentinel-1获取的多期数据，以更高分辨率的GF-3影像为提取对象，以同期Sentinel-2影像作为辅助进行校正，应用成熟的阈值分割法，提取了多期洪水淹没范围(图6)。

图6 研究区域SAR卫星多时相水体提取结果

同时，根据不同覆盖类型对未受灾年多时相NDVI数据进行非监督分类。通过地面调查和获取到的受灾年5月28日Sentinel2光学影像(10m分辨率)的目视解译，结合土地覆盖分类体系不同地类NDVI多时相变化的季节性、峰值特征，分为农田、森林、草地、灌木丛、湿地、水体、建筑用地和荒地8个类型。将上述受灾年提取的水体信息与聚类结果叠加图(图7)同地类被淹没的面积随时相的变化，由于湿地和荒地所占比例过小忽略不计。

图7 研究区域聚类结果与各地类动态信息统计图

(一)农田淹没面积对比验证

由GF-3和Sentinel1卫星数据提取到的5月30日农田被淹没面积为22.05km2,MODIS则取对应于5月30日GF-3和Sentinel1时相提取的2017年与2016年的NDVI差值从该时相开始连续3个时相小大于0.18时提取的像元数。计算得出淹没面积为25.19km2略大于GF-3卫星提取结果，其原因可能是NDVI时序数据提取的特征像元有来自已经退水但NDVI值无法恢复到未受灾年的区域。同时由于MODIS影像分辨率较低，无法精确提取微小水面信息，从而平滑了那些小面积的水体变化信息。

(二)农田淹没范围

由GF-3和Sentinel1灾中和灾后多时相数据计算出农田被淹没的面积随时相的变化可知农田部分至8月8日水淹面积为4.23km2，相较于灾中已经褪去82.6%。在MODIS时序中统计2017受灾年持续时间Dn小于65天既持续影响时间在8月8日前后的面积为19.7km2，占整个持续天数比例的78.6%，两组数据基本一致。MODIS占比相对较少的原因可能是在水退之后由于土壤过饱和会继续影响作物生长。可以反映出大部分区域受灾害影响的持续时间在两个月左右。

结束语

洪灾通过使农田土壤水分过饱和从而影响农作物的生长是一个渐变的过程。通过多源遥感卫星的多时相的数据进行动态监测才能准确评估其受灾程度和灾后恢复情况。利用获取到的Terra/MODIS受灾年和未受灾年观测获取的多时相NDVI数据的变化特征，提出了通过设定阈值来提取灾害的动态信息的方法。并利用GF-3和Sentinel1卫星对受灾区域进行辅助，实现了对该区域更精确的动态监测。结果表明，基于MODIS多时相合成的NDVI时序数据可提取不同受灾区域的开始时间、持续时间以及灾后恢复状况等诸多时间和空间的动态信息。通过高分辨率微波遥感卫星提取的水体信息进行对比验证后发现该方法可以快速准确的获取大范围灾情信息。为研究暴雨、洪灾以及其他覆盖面积较大的突发性地质灾害提供了一个方法应用研究案例。

通讯作者：张文君