基于多源衛星的洪水災害動態信息提取

周 帆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時序數據可提取不同受災區域的開始時間、持續時間以及災后恢復狀況等諸多時間和空間的動態信息。通過高分辨率微波遙感衛星提取的水體信息進行對比驗證后發現該方法可以快速準確的獲取大范圍災情信息。為研究暴雨、洪災以及其他覆蓋面積較大的突發性地質災害提供了一個方法應用研究案例。

通訊作者：張文君