耕地暴雨灾情动态信息的遥感提取

郭旭东，苏亚丽，吴长江，陈瑜琦，雷莉萍，廖静娟，汪晓帆，张衍毓

（1.中国土地勘测规划院 国土资源部土地利用重点实验室，北京 100035；2.中国科学院遥感与数字地球研究所数字地球重点实验室，北京 100094；3.西安科技大学 研究生院，陕西 西安 710054）

在全球气候变暖背景下，持续降水以及暴雨等气象灾害频繁发生，严重影响了农业生产活动。暴雨灾害是指在地势低洼、地形闭塞的地区，一次短时的或连续的强降水过程中雨水不能迅速宣泄造成农田积水和土壤水分过度饱和而给农业带来的灾害。由于暴雨发生时间和强降水持续时长的不同，导致耕地的受灾程度、面积以及后期农作物恢复情况均具有随时间变化的动态特征。这些特征是农业损失评估和国家灾害救助所需的重要灾情信息。

能够客观准确地获取地面真实状况以及刻画地表动态过程的卫星遥感观测技术已成为灾情评估的重要手段[1]。目前多数研究关注于利用卫星遥感观测数据提取暴雨洪水导致的耕地淹没空间信息[2-6]，而现实中暴雨和持续强降水过程所引起的灾害，不仅包含耕地淹没范围，还包含从耕地受灾发生的时间、受害程度到灾后恢复的动态过程，因此不仅在空间上而且在时间上也需获取灾情的动态信息。

本文利用卫星遥感观测的时序数据，结合对应暴雨灾害发生的时间、范围、受害程度以及灾后恢复状况的时序变化特征，研发了受灾信息提取方法，为作物灾后损失评估和国家灾害救助提供了重要的决策依据。

1 数据来源与研究方法

1.1 NDVI时序数据

本文以安徽省巢湖地区为实验区，收集了Terra/MODIS观测获取的NDVI时间序列数据。2016 年夏季巢湖地区发生的持续强降水和暴雨天气导致大面积耕地被淹和农作物受灾，因此本文收集了2016年和2015 年的MODIS观测数据处理的16 d合成NDVI数据产品。按16 d的周期每年有23个时相，但由于年末最后一个时相缺失，2016年时相数为22个。该数据产品来源于美国NASA的EOS系列卫星Terra观测获取的植被指数产品，全称为MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid (MOD13Q1)[7]。本文收集的NDVI数据产品是该植被指数产品之一，空间分辨率为250 m。16 d合成的NDVI数据是由最佳像元选择合成算法处理得到的数据产品[8-9]。该算法以16 d为单位，在一定传感器观测角范围内，从16 d中选出无云或受云影响最小（即NDVI值最大）、且距星下点最近的像元作为最佳像元，通过逐像元选择处理16 d的最佳像元，合成得到一幅NDVI数据。

1.2 NDVI时序变化的灾害特征

持续强降水导致耕地土壤水分过度饱和，暴雨甚至淹没耕地，必将影响耕地作物的生长，此时卫星观测的地表反射率也必将随耕地土壤和作物的变化而发生变化，表现为NDVI值的持续降低，因此可利用NDVI变化与耕地受灾后变化的关系特征提取灾情信息。

以巢湖地区2016年暴雨受灾区域为例，结合地面调查，选择6个不同受灾程度样本区域，提取了2016 年和未受灾的2015年的NDVI多时相数据（图1），并对比其变化特征。

1）不同样本区受灾年的NDVI值开始降低的时期不同。对比研究区合肥市气象观测站2016年的降水量变化（图2）可知，在暴雨发生的7月1日之前，2016 年的NDVI值已开始下降；图1a、1b样本区已于5月底、6月初开始下降；图1c～1e样本区于6月底开始下降，表明在暴雨发生前的连续降水已导致土壤水分饱和，并开始影响农作物的生长。

图1 2016年(红色)和2015年(蓝色)的NDVI多时相变化特征

2）不同样本区受灾后NDVI值降低的持续时间不同。2016年开始受灾时期的NDVI值比2015年同时期减小约0.15，以此为阈值统计2015年与2016年同时相的NDVI差值连续大于0.15的时相个数，图 1a～1c分别呈现了9、6、8期的时相个数，这些样本区的NDVI值在受灾时降低，随后又回升到了2015年的NDVI值。图1d、1e样本区分别连续有11、12个时相的NDVI差值大于0.15，其NDVI值受灾降低后直至年末也没有回升，这与2016年7月1日暴雨后的连续降水导致后期没有进行作物耕种有关。低NDVI值持续的时相个数揭示了耕地受灾后作物耕种生长恢复的情况，可作为一 个特征值评估耕地恢复的状态，即受灾程度。

3）非灾害所致NDVI值的降低。2016年NDVI值的减小除灾害原因外，也还可能有其他变化的因素，如耕地变为建筑用地、休耕地等。如图1f所示，2016 年NDVI值从起始点开始就显示了比2015年低的特征，且没有季节性变化，经地面调查核实发现，该区域的土地利用在2016年由耕地改为了水产养殖田。

图2 合肥市气象站2016年降水量变化

1.3 灾情动态信息提取方法

根据对持续强降水和暴雨导致的NDVI值减小的特征分析，本文提出了利用NDVI时序变化提取耕地受灾动态信息方法（图3）。其主要流程包括异常数据的剔除及其插补、NDVI时序变化检测、灾害动态信息提取、受灾时期和时长识别以及灾后恢复状况分类。

图3 基于NDVI时序数据的暴雨灾害耕地受灾动态信息提取方法流程图

1）异常数据的剔除及其插补。MODIS 16 d NDVI产品数据虽由最佳像元合成方法处理得到，很大程度上减小了云的影响[8]，但由于中国南方地区天气多云雨，会出现16 d均为多云天气的情况，导致NDVI合成数据产品中仍存在一些受云影响的数据，使得NDVI的多时相变化异于正常的变化规律。

在NDVI多时相变化中，农作物从生长到成熟收割，NDVI呈先增大后减小的变化规律[9]。即使NDVI值在收割或耕地受灾后会减小，但不会在后一时相出现随之增大的变化特征，因此当NDVI值与前后相邻时相比较时，若同时呈突然减小或增大的变化特征，则可推测其为受云等影响的异常值。根据NDVI值的正常多时相变化特征，当像元NDVI值与前后两个时相的绝对差值均大于0.15时，判断其为异常值，并利用前后时相NDVI值进行线性插值，由此减小异常值对多时相变化特征分析的影响。

2）变化检测与灾情动态信息的提取。利用受灾年和未受灾年NDVI的差值，根据土地利用辅助数据提取耕地像元，再逐像元进行多时相的统计分析提取受灾动态信息。

令像元每年NDVI的时相总数为N，计算时相（i）未受灾年（NDVIc[i]）与受灾年（NDVIs

[i]）的差值ΔNDVI[i]。设定由于暴雨灾害导致地表积水与土壤水分增加引起的NDVI降低的差值阈值为T，若ΔNDVI[i]连续2个时相以上均大于T，则以时相i为受灾开始时相（Dt），累计大于阈值的时相个数，一旦小于阈值停止累计，输出满足条件的时相个数（Dn），该时相个数可表征地表受灾害影响后持续的时长。

根据我国只有在春季4月之后（i=6）才有可能发生持续强降水和暴雨的规律，以Dt必大于6或Dn必小于（N-6）为条件限制，可剔除图1f中由非灾害所引起的NDVI值变化特征。

2 研究结果与讨论

2.1 应用处理结果

本文利用巢湖地区2016年和未受灾的2015年MODIS多时相NDVI数据对上述方法进行实验验证。首先分别对实验区收集的2 a MODIS数据进行异常剔除和插补处理，统一得到21个时相的NDVI数据；再根据2016年和2015年多时相NDVI的变化检测，提取巢湖地区耕地暴雨灾情动态信息。

图4的开始时间是根据NDVI数据的每个时相间隔16 d，将开始时相Dt进行日期换算后的结果。图5的持续天数也是根据每个时相间隔16 d，由时相个数Dn换算的天数。由受灾提取结果可知，在整个研究区（左上图），耕地受灾主要发生在巢湖东南部的低洼区域。大部分耕地受灾开始在暴雨发生的6月25日～7 月26日（图4中黄色和青色区域）；但在强大暴雨发生的7月之前，由于持续强降水引起的土壤水分饱和已开始影响了农作物的生长，导致2016年的一些耕地的NDVI值比2015年小，受灾开始时间显示在强暴雨之前的6月（图4中红色区域），说明了该区域低洼容易积水受灾的特征。图5的暴雨受灾持续天数主要显示在48～112 d，持续天数越多说明受灾害影响越大。以6月开始受灾计算，显示在160 d和176 d的结果表明，该区域受灾后其NDVI没有回到2015年的水平，说明受灾后基本没有恢复。

研究区受灾总面积为1 193.2 km2，其中近80%分布于巢湖东南地区的巢湖市、庐江县、桐城市、枞阳县和无为县；受灾区域的55%（面积为665 km2）显示为受灾持续天数少于48 d（图5中绿色和青色区域），表明该区域耕地受害程度轻；27%的区域（面积为324 km2）受灾持续天数约为112 d（图5中黄色），即约3个多月地表状态没有恢复到未受灾年的状态；18%的区域（面积为223 km2）受灾持续天数为144 d以上（图 5中红色和紫红区域），即以6月开始受灾计算，该区域到11月还未恢复到2015年的状态，表明受灾严重，耕地作物的种植生长未恢复到2015年的状态。

图4 2016年受灾开始时间的识别结果

图5 2016年受灾持续天数的计算结果

2.2 结果对比与验证

1）研究区NDVI多时相变化特征的聚类。利用NDVI的多时相变化特征，可以很好地获得耕地农作物类型[10]。根据林地、农作物、水体、居民地等的NDVI季节变化特征，对1 a的多时相NDVI数据进行非监督分类。图6为2016年和2015年非监督分类结果以及对应类型的各时相NDVI均值的多时相变化特征。根据地面调查和参考卫星影像的目视判别，结合各类NDVI值及其多时相变化的季节性、峰值等特征，将聚类结果归纳为5个类型，分别为耕地-1 （1季作物）、耕地-2（2季作物）、林地、水体和其他/居民地，这里的其他为NDVI值较低且季节变化幅度较小的居民地，如撂荒地、草地、河滩等。

图6 NDVI季节变化特征的非监督聚类结果及其对应的各类NDVI均值的时序变化

对比2016年和2015年的分类结果可知，两年巢湖北部的各类型空间分布非常相似，但在巢湖的东南部区域（图6a、6b中圆内区域）可以发现，2015年很多耕地（图6中绿色和黄色区域）在2016年变化为其他/居民地（图6中桔黄色区域）和水体（图6中蓝色区域）。研究区范围内各类型面积统计结果显示，与2015年相比，2016年的林地（图6中青色区域）面积没有变化，但水体和其他的面积分别增大了26%和12%，耕地减小了15%，这与2016年的持续强降水和暴雨天气导致的耕地受淹有关。对比图4、5可知，这些变化区域与受灾区域的分布一致。

2）与其他卫星遥感数据提取结果的对比。收集研究区暴雨发生后的2016年7月24日Sentinel-1观测获取的雷达数据（空间分辨率为10 m），利用常规方法提取水淹区域，并与本文方法提取的结果进行比较。在空间分布上，图4、5显示的基于NDVI时序变化提取的受灾区域与Sentinel-1的结果一致。Sentinel-1区域的受淹面积为1 082 km2（图7），对应区域图4中开始时期8月11日以前的受灾面积为1 178 km2，略大于Sentinel-1结果，相差9%。

图7 Sentinel-1雷达数据提取的受灾区域（红色）图

3 结 语

卫星遥感观测已经能提供大量且具有一定精度的多年时间序列植被指数数据，这些数据可揭示农作物生长的季节变化特征和耕地状态信息。暴雨天气以及持续强降水导致的土壤水分饱和，对农作物生长的影响是一个动态过程，因此本文利用受灾和未受灾的MODIS-NDVI多时序数据，提出了一种计算方便、处理简单的灾情信息遥感提取方法。该方法不同于仅提取暴雨洪水淹没空间信息的常规方法，不但可提取大范围区域的受灾空间信息而且能获取由持续强降水和暴雨天气灾害导致的耕地受灾时间的动态信息以及灾后恢复状况信息。这些灾情动态信息能更好地为灾后耕地的恢复以及国家灾后损失评估和救助决策提供科学依据。

[1] HU Z W, GONG H L, ZHU L Y. Fast Flooding Information Extraction in Emergency Response of Flood Disaster[J].International Society for Photogrammetry and Remote Sensing,2012,54(4):173-177

[2] 许超,蒋卫国,万立冬,等.基于MODIS时间序列数据的洞庭湖区洪水淹没频率研究[J].灾害学,2016,31(1):96-101

[3] Shrestha R, DI L P, YU G N, et al. Detection of Flood and Its Impact on Crops Using NDVI-corn Case[C].Second International Conference on Agro-geoinformatics,2013,7(2):200-204

[4] Kuenzer C, GUO H D, Huth J, et al. Flood Mapping and Flood Dynamics of the Mekong Delta: ENVISAT-ASAR-WSM Based Time Series Analyses[J]. Remote Sensing,2013,5(2):687-715

[5] Sakamoto T, Nguyen N V, Kotera A, et al. Detecting Temporal Changes in the Extent of Annual Flooding Within the Cambodia and the Vietnamese Mekong Delta from MODIS Time-series Imagery[J]. Remote Sensing of Environment,2007,109(3):295-313

[6] 李通,张丽,申茜,等.湄公河下游洪灾淹没面积多源遥感时序监测分析[J].应用科学学报,2016,34(1):75-83

[7] 美国国家航空航天局(NASA)陆地过程分布式数据档案中心.MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid[EB/OL].[2016-12-02]. http://ladsweb.nascom.nasa.gov/data/search.html

[8] National Aeronautics and Space Administration.MODIS Vegetation Index User's Guide(MOD13 Series)[EB/OL].[2015-06-03]. http://modis-land.gsfc.nasa.gov/vi.html

[9] Leeuwen W J D V, Huete A R, Laing T W. MODIS Vegetation Index Compositing Approach: a Prototype with AVHRR Data [J].Remote Sensing of Environment,1999,69(3):264-280

[10] JIANG D, WANG N B, YANG X H, et al. Study on the Interaction Between NDVI Profile and the Growing Status of Crops[J]. Chinese Geographical Science,2003,13(1):62-65

[11] 梁益同,刘可群,周守华,等.EOS-MODIS数据监测暴雨洪涝灾害的技术方法[J].暴雨灾害,2008,27(1):64-67

[12] 张焕雪,曹欣,李强子,等.基于多时相环境星NDVI时间序列的农作物分类研究[J].遥感技术与应用,2015,30(2):304-311

[13] 刘寒,冯莉,朱榴骏,等.环境星归一化植被指数时间序列滤波算法比较[J].遥感信息,2015,30(5):69-76