郭旭東，蘇亞麗，吳長江，陳瑜琦，雷莉萍，廖靜娟，汪曉帆，張衍毓

（1.中國土地勘測規(guī)劃院 國土資源部土地利用重點實驗室，北京 100035；2.中國科學院遙感與數(shù)字地球研究所數(shù)字地球重點實驗室，北京 100094；3.西安科技大學 研究生院，陜西 西安 710054）

在全球氣候變暖背景下，持續(xù)降水以及暴雨等氣象災害頻繁發(fā)生，嚴重影響了農(nóng)業(yè)生產(chǎn)活動。暴雨災害是指在地勢低洼、地形閉塞的地區(qū)，一次短時的或連續(xù)的強降水過程中雨水不能迅速宣泄造成農(nóng)田積水和土壤水分過度飽和而給農(nóng)業(yè)帶來的災害。由于暴雨發(fā)生時間和強降水持續(xù)時長的不同，導致耕地的受災程度、面積以及后期農(nóng)作物恢復情況均具有隨時間變化的動態(tài)特征。這些特征是農(nóng)業(yè)損失評估和國家災害救助所需的重要災情信息。

能夠客觀準確地獲取地面真實狀況以及刻畫地表動態(tài)過程的衛(wèi)星遙感觀測技術(shù)已成為災情評估的重要手段[1]。目前多數(shù)研究關(guān)注于利用衛(wèi)星遙感觀測數(shù)據(jù)提取暴雨洪水導致的耕地淹沒空間信息[2-6]，而現(xiàn)實中暴雨和持續(xù)強降水過程所引起的災害，不僅包含耕地淹沒范圍，還包含從耕地受災發(fā)生的時間、受害程度到災后恢復的動態(tài)過程，因此不僅在空間上而且在時間上也需獲取災情的動態(tài)信息。

本文利用衛(wèi)星遙感觀測的時序數(shù)據(jù)，結(jié)合對應暴雨災害發(fā)生的時間、范圍、受害程度以及災后恢復狀況的時序變化特征，研發(fā)了受災信息提取方法，為作物災后損失評估和國家災害救助提供了重要的決策依據(jù)。

1 數(shù)據(jù)來源與研究方法

1.1 NDVI時序數(shù)據(jù)

本文以安徽省巢湖地區(qū)為實驗區(qū)，收集了Terra/MODIS觀測獲取的NDVI時間序列數(shù)據(jù)。2016 年夏季巢湖地區(qū)發(fā)生的持續(xù)強降水和暴雨天氣導致大面積耕地被淹和農(nóng)作物受災，因此本文收集了2016年和2015 年的MODIS觀測數(shù)據(jù)處理的16 d合成NDVI數(shù)據(jù)產(chǎn)品。按16 d的周期每年有23個時相，但由于年末最后一個時相缺失，2016年時相數(shù)為22個。該數(shù)據(jù)產(chǎn)品來源于美國NASA的EOS系列衛(wèi)星Terra觀測獲取的植被指數(shù)產(chǎn)品，全稱為MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid (MOD13Q1)[7]。本文收集的NDVI數(shù)據(jù)產(chǎn)品是該植被指數(shù)產(chǎn)品之一，空間分辨率為250 m。16 d合成的NDVI數(shù)據(jù)是由最佳像元選擇合成算法處理得到的數(shù)據(jù)產(chǎn)品[8-9]。該算法以16 d為單位，在一定傳感器觀測角范圍內(nèi)，從16 d中選出無云或受云影響最?。碞DVI值最大）、且距星下點最近的像元作為最佳像元，通過逐像元選擇處理16 d的最佳像元，合成得到一幅NDVI數(shù)據(jù)。

1.2 NDVI時序變化的災害特征

持續(xù)強降水導致耕地土壤水分過度飽和，暴雨甚至淹沒耕地，必將影響耕地作物的生長，此時衛(wèi)星觀測的地表反射率也必將隨耕地土壤和作物的變化而發(fā)生變化，表現(xiàn)為NDVI值的持續(xù)降低，因此可利用NDVI變化與耕地受災后變化的關(guān)系特征提取災情信息。

以巢湖地區(qū)2016年暴雨受災區(qū)域為例，結(jié)合地面調(diào)查，選擇6個不同受災程度樣本區(qū)域，提取了2016 年和未受災的2015年的NDVI多時相數(shù)據(jù)（圖1），并對比其變化特征。

1）不同樣本區(qū)受災年的NDVI值開始降低的時期不同。對比研究區(qū)合肥市氣象觀測站2016年的降水量變化（圖2）可知，在暴雨發(fā)生的7月1日之前，2016 年的NDVI值已開始下降；圖1a、1b樣本區(qū)已于5月底、6月初開始下降；圖1c～1e樣本區(qū)于6月底開始下降，表明在暴雨發(fā)生前的連續(xù)降水已導致土壤水分飽和，并開始影響農(nóng)作物的生長。

圖1 2016年(紅色)和2015年(藍色)的NDVI多時相變化特征

2）不同樣本區(qū)受災后NDVI值降低的持續(xù)時間不同。2016年開始受災時期的NDVI值比2015年同時期減小約0.15，以此為閾值統(tǒng)計2015年與2016年同時相的NDVI差值連續(xù)大于0.15的時相個數(shù)，圖 1a～1c分別呈現(xiàn)了9、6、8期的時相個數(shù)，這些樣本區(qū)的NDVI值在受災時降低，隨后又回升到了2015年的NDVI值。圖1d、1e樣本區(qū)分別連續(xù)有11、12個時相的NDVI差值大于0.15，其NDVI值受災降低后直至年末也沒有回升，這與2016年7月1日暴雨后的連續(xù)降水導致后期沒有進行作物耕種有關(guān)。低NDVI值持續(xù)的時相個數(shù)揭示了耕地受災后作物耕種生長恢復的情況，可作為一 個特征值評估耕地恢復的狀態(tài)，即受災程度。

3）非災害所致NDVI值的降低。2016年NDVI值的減小除災害原因外，也還可能有其他變化的因素，如耕地變?yōu)榻ㄖ玫?、休耕地等。如圖1f所示，2016 年NDVI值從起始點開始就顯示了比2015年低的特征，且沒有季節(jié)性變化，經(jīng)地面調(diào)查核實發(fā)現(xiàn)，該區(qū)域的土地利用在2016年由耕地改為了水產(chǎn)養(yǎng)殖田。

圖2 合肥市氣象站2016年降水量變化

1.3 災情動態(tài)信息提取方法

根據(jù)對持續(xù)強降水和暴雨導致的NDVI值減小的特征分析，本文提出了利用NDVI時序變化提取耕地受災動態(tài)信息方法（圖3）。其主要流程包括異常數(shù)據(jù)的剔除及其插補、NDVI時序變化檢測、災害動態(tài)信息提取、受災時期和時長識別以及災后恢復狀況分類。

圖3 基于NDVI時序數(shù)據(jù)的暴雨災害耕地受災動態(tài)信息提取方法流程圖

1）異常數(shù)據(jù)的剔除及其插補。MODIS 16 d NDVI產(chǎn)品數(shù)據(jù)雖由最佳像元合成方法處理得到，很大程度上減小了云的影響[8]，但由于中國南方地區(qū)天氣多云雨，會出現(xiàn)16 d均為多云天氣的情況，導致NDVI合成數(shù)據(jù)產(chǎn)品中仍存在一些受云影響的數(shù)據(jù)，使得NDVI的多時相變化異于正常的變化規(guī)律。

在NDVI多時相變化中，農(nóng)作物從生長到成熟收割，NDVI呈先增大后減小的變化規(guī)律[9]。即使NDVI值在收割或耕地受災后會減小，但不會在后一時相出現(xiàn)隨之增大的變化特征，因此當NDVI值與前后相鄰時相比較時，若同時呈突然減小或增大的變化特征，則可推測其為受云等影響的異常值。根據(jù)NDVI值的正常多時相變化特征，當像元NDVI值與前后兩個時相的絕對差值均大于0.15時，判斷其為異常值，并利用前后時相NDVI值進行線性插值，由此減小異常值對多時相變化特征分析的影響。

2）變化檢測與災情動態(tài)信息的提取。利用受災年和未受災年NDVI的差值，根據(jù)土地利用輔助數(shù)據(jù)提取耕地像元，再逐像元進行多時相的統(tǒng)計分析提取受災動態(tài)信息。

令像元每年NDVI的時相總數(shù)為N，計算時相（i）未受災年（NDVIc[i]）與受災年（NDVIs

[i]）的差值ΔNDVI[i]。設(shè)定由于暴雨災害導致地表積水與土壤水分增加引起的NDVI降低的差值閾值為T，若ΔNDVI[i]連續(xù)2個時相以上均大于T，則以時相i為受災開始時相（Dt），累計大于閾值的時相個數(shù)，一旦小于閾值停止累計，輸出滿足條件的時相個數(shù)（Dn），該時相個數(shù)可表征地表受災害影響后持續(xù)的時長。

根據(jù)我國只有在春季4月之后（i=6）才有可能發(fā)生持續(xù)強降水和暴雨的規(guī)律，以Dt必大于6或Dn必小于（N-6）為條件限制，可剔除圖1f中由非災害所引起的NDVI值變化特征。

2 研究結(jié)果與討論

2.1 應用處理結(jié)果

本文利用巢湖地區(qū)2016年和未受災的2015年MODIS多時相NDVI數(shù)據(jù)對上述方法進行實驗驗證。首先分別對實驗區(qū)收集的2 a MODIS數(shù)據(jù)進行異常剔除和插補處理，統(tǒng)一得到21個時相的NDVI數(shù)據(jù)；再根據(jù)2016年和2015年多時相NDVI的變化檢測，提取巢湖地區(qū)耕地暴雨災情動態(tài)信息。

圖4的開始時間是根據(jù)NDVI數(shù)據(jù)的每個時相間隔16 d，將開始時相Dt進行日期換算后的結(jié)果。圖5的持續(xù)天數(shù)也是根據(jù)每個時相間隔16 d，由時相個數(shù)Dn換算的天數(shù)。由受災提取結(jié)果可知，在整個研究區(qū)（左上圖），耕地受災主要發(fā)生在巢湖東南部的低洼區(qū)域。大部分耕地受災開始在暴雨發(fā)生的6月25日～7 月26日（圖4中黃色和青色區(qū)域）；但在強大暴雨發(fā)生的7月之前，由于持續(xù)強降水引起的土壤水分飽和已開始影響了農(nóng)作物的生長，導致2016年的一些耕地的NDVI值比2015年小，受災開始時間顯示在強暴雨之前的6月（圖4中紅色區(qū)域），說明了該區(qū)域低洼容易積水受災的特征。圖5的暴雨受災持續(xù)天數(shù)主要顯示在48～112 d，持續(xù)天數(shù)越多說明受災害影響越大。以6月開始受災計算，顯示在160 d和176 d的結(jié)果表明，該區(qū)域受災后其NDVI沒有回到2015年的水平，說明受災后基本沒有恢復。

研究區(qū)受災總面積為1 193.2 km2，其中近80%分布于巢湖東南地區(qū)的巢湖市、廬江縣、桐城市、樅陽縣和無為縣；受災區(qū)域的55%（面積為665 km2）顯示為受災持續(xù)天數(shù)少于48 d（圖5中綠色和青色區(qū)域），表明該區(qū)域耕地受害程度輕；27%的區(qū)域（面積為324 km2）受災持續(xù)天數(shù)約為112 d（圖5中黃色），即約3個多月地表狀態(tài)沒有恢復到未受災年的狀態(tài)；18%的區(qū)域（面積為223 km2）受災持續(xù)天數(shù)為144 d以上（圖 5中紅色和紫紅區(qū)域），即以6月開始受災計算，該區(qū)域到11月還未恢復到2015年的狀態(tài)，表明受災嚴重，耕地作物的種植生長未恢復到2015年的狀態(tài)。

圖4 2016年受災開始時間的識別結(jié)果

圖5 2016年受災持續(xù)天數(shù)的計算結(jié)果

2.2 結(jié)果對比與驗證

1）研究區(qū)NDVI多時相變化特征的聚類。利用NDVI的多時相變化特征，可以很好地獲得耕地農(nóng)作物類型[10]。根據(jù)林地、農(nóng)作物、水體、居民地等的NDVI季節(jié)變化特征，對1 a的多時相NDVI數(shù)據(jù)進行非監(jiān)督分類。圖6為2016年和2015年非監(jiān)督分類結(jié)果以及對應類型的各時相NDVI均值的多時相變化特征。根據(jù)地面調(diào)查和參考衛(wèi)星影像的目視判別，結(jié)合各類NDVI值及其多時相變化的季節(jié)性、峰值等特征，將聚類結(jié)果歸納為5個類型，分別為耕地-1 （1季作物）、耕地-2（2季作物）、林地、水體和其他/居民地，這里的其他為NDVI值較低且季節(jié)變化幅度較小的居民地，如撂荒地、草地、河灘等。

圖6 NDVI季節(jié)變化特征的非監(jiān)督聚類結(jié)果及其對應的各類NDVI均值的時序變化

對比2016年和2015年的分類結(jié)果可知，兩年巢湖北部的各類型空間分布非常相似，但在巢湖的東南部區(qū)域（圖6a、6b中圓內(nèi)區(qū)域）可以發(fā)現(xiàn)，2015年很多耕地（圖6中綠色和黃色區(qū)域）在2016年變化為其他/居民地（圖6中桔黃色區(qū)域）和水體（圖6中藍色區(qū)域）。研究區(qū)范圍內(nèi)各類型面積統(tǒng)計結(jié)果顯示，與2015年相比，2016年的林地（圖6中青色區(qū)域）面積沒有變化，但水體和其他的面積分別增大了26%和12%，耕地減小了15%，這與2016年的持續(xù)強降水和暴雨天氣導致的耕地受淹有關(guān)。對比圖4、5可知，這些變化區(qū)域與受災區(qū)域的分布一致。

2）與其他衛(wèi)星遙感數(shù)據(jù)提取結(jié)果的對比。收集研究區(qū)暴雨發(fā)生后的2016年7月24日Sentinel-1觀測獲取的雷達數(shù)據(jù)（空間分辨率為10 m），利用常規(guī)方法提取水淹區(qū)域，并與本文方法提取的結(jié)果進行比較。在空間分布上，圖4、5顯示的基于NDVI時序變化提取的受災區(qū)域與Sentinel-1的結(jié)果一致。Sentinel-1區(qū)域的受淹面積為1 082 km2（圖7），對應區(qū)域圖4中開始時期8月11日以前的受災面積為1 178 km2，略大于Sentinel-1結(jié)果，相差9%。

圖7 Sentinel-1雷達數(shù)據(jù)提取的受災區(qū)域（紅色）圖

3 結(jié) 語

衛(wèi)星遙感觀測已經(jīng)能提供大量且具有一定精度的多年時間序列植被指數(shù)數(shù)據(jù)，這些數(shù)據(jù)可揭示農(nóng)作物生長的季節(jié)變化特征和耕地狀態(tài)信息。暴雨天氣以及持續(xù)強降水導致的土壤水分飽和，對農(nóng)作物生長的影響是一個動態(tài)過程，因此本文利用受災和未受災的MODIS-NDVI多時序數(shù)據(jù)，提出了一種計算方便、處理簡單的災情信息遙感提取方法。該方法不同于僅提取暴雨洪水淹沒空間信息的常規(guī)方法，不但可提取大范圍區(qū)域的受災空間信息而且能獲取由持續(xù)強降水和暴雨天氣災害導致的耕地受災時間的動態(tài)信息以及災后恢復狀況信息。這些災情動態(tài)信息能更好地為災后耕地的恢復以及國家災后損失評估和救助決策提供科學依據(jù)。

[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] 許超,蔣衛(wèi)國,萬立冬,等.基于MODIS時間序列數(shù)據(jù)的洞庭湖區(qū)洪水淹沒頻率研究[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] 李通,張麗,申茜,等.湄公河下游洪災淹沒面積多源遙感時序監(jiān)測分析[J].應用科學學報,2016,34(1):75-83

[7] 美國國家航空航天局(NASA)陸地過程分布式數(shù)據(jù)檔案中心.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數(shù)據(jù)監(jiān)測暴雨洪澇災害的技術(shù)方法[J].暴雨災害,2008,27(1):64-67

[12] 張煥雪,曹欣,李強子,等.基于多時相環(huán)境星NDVI時間序列的農(nóng)作物分類研究[J].遙感技術(shù)與應用,2015,30(2):304-311

[13] 劉寒,馮莉,朱榴駿,等.環(huán)境星歸一化植被指數(shù)時間序列濾波算法比較[J].遙感信息,2015,30(5):69-76