基于多源衛(wèi)星的洪水災(zāi)害動態(tài)信息提取

周 帆1,2張文君1 雷莉萍2 張智杰2 郭曉東 汪曉帆

(1.西南科技大學(xué)環(huán)境與資源學(xué)院 四川 綿陽 621010)(2.中國科學(xué)院遙感與數(shù)字地球研究所數(shù)字地球重點實驗室 北京 100094)(3.中國土地勘測規(guī)劃院國土資源部土地利用重點實驗室 北京 100035)

引言

目前，洪水災(zāi)害被認(rèn)為是世界上危害性最大的自然災(zāi)害之一，其造成的損失在全球所有自然災(zāi)害中占40%。洪水往往發(fā)生在人口密度高，河流和湖泊集中，降雨量大，農(nóng)業(yè)種植的程度高的地方[1]。其對農(nóng)業(yè)的影響主要表現(xiàn)為淹沒耕地的大量的積水短時間內(nèi)難以褪去，農(nóng)田的受災(zāi)面積和災(zāi)后恢復(fù)狀況均隨著洪災(zāi)發(fā)生和褪去的時間不同而反應(yīng)出不同的動態(tài)特征。這些特征是對災(zāi)情動態(tài)監(jiān)測以及災(zāi)后損失評估的重要依據(jù)。

由于洪災(zāi)常伴隨著持續(xù)的強(qiáng)降雨，具有突發(fā)性而不易實地調(diào)查。傳統(tǒng)的調(diào)查和監(jiān)測方法無法取得實時可靠的數(shù)據(jù)。遙感技術(shù)因其大范圍、全天候全天時的特點被廣泛應(yīng)用于洪災(zāi)的監(jiān)測，利用多時相的遙感影像可以實時監(jiān)測洪災(zāi)的發(fā)生狀況。喻光明等人分析了四湖流域“91.7”暴雨分區(qū)水量平衡及地面水情實況,利用遙感技術(shù)研究了該區(qū)“91.7”暴雨的洪水勢態(tài),包括洪水淹沒范圍估算及其時空特性分析,據(jù)此作出了災(zāi)情評價[2]。張國平等人利用遙感,對多種地表水體如河流和湖泊進(jìn)行空間識別、定位及定量計算面積等研究內(nèi)容開展洪水水情監(jiān)測[3]。王云秀等人利用高分與資源衛(wèi)星對暴雨后河北省中南部水庫面積進(jìn)行了監(jiān)測分析[4]。李健等人利用利用FY-3A、HJ-1A/B和EOS多源衛(wèi)星遙感數(shù)據(jù),結(jié)合地面氣象觀測數(shù)據(jù)和基礎(chǔ)地理信息數(shù)據(jù),對2010年7月下旬至8月初吉林省特大暴雨導(dǎo)致的洪澇災(zāi)害進(jìn)行了監(jiān)測[5]。但在現(xiàn)實中洪災(zāi)對區(qū)域產(chǎn)生的影響不僅有農(nóng)田淹沒范圍信息也應(yīng)分析其受災(zāi)發(fā)生的時間、持續(xù)時間以及災(zāi)后恢復(fù)的動態(tài)信息。

本文利用多源遙感衛(wèi)星的多時相數(shù)據(jù)對洪災(zāi)發(fā)生的時間、范圍、不同區(qū)域的受災(zāi)程度以及災(zāi)后恢復(fù)狀況進(jìn)行動態(tài)監(jiān)測。通過時序變化檢測提取農(nóng)田受災(zāi)動態(tài)信息的方法，為受災(zāi)狀況分析以及災(zāi)后恢復(fù)評估提供決策的重要科學(xué)依據(jù)。

一、數(shù)據(jù)采集

由于2017年5月25日開始的季風(fēng)性降雨導(dǎo)致了斯里蘭卡境內(nèi)山洪暴發(fā)，這次水災(zāi)是自2004年印度洋海嘯以來，斯里蘭卡遭遇的最大一次自然災(zāi)害。我們以洪水暴發(fā)的集中區(qū)域斯里蘭卡南部馬塔勒區(qū)作為研究區(qū)域(圖1)。

圖1 研究區(qū)域的地理位置和衛(wèi)星遙感彩色影像

(2016年1月29日Landsat-8衛(wèi)星觀測的可見光和紅外光(OLI6-5-3波段)波段合成的影像)

本文收集了2017年和未受災(zāi)的前一年2016年MODIS觀測數(shù)據(jù)處理的16天合成NDVI數(shù)據(jù)產(chǎn)品(表1)。由于研究區(qū)域處于不同的軌道，將數(shù)據(jù)拼接后以16天為周期的范圍內(nèi)由受云影響最小以及距離星下最近的最佳像元合成，每年有23個時相。該數(shù)據(jù)產(chǎn)品來源于美國NASA的EOS系列衛(wèi)星Terra觀測獲取的植被指數(shù)產(chǎn)品，全稱為MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid(MOD13Q1)[6]。本文收集的NDVI數(shù)據(jù)產(chǎn)品是該植被指數(shù)產(chǎn)品之一，數(shù)據(jù)的空間分辨率為250m。16天合成的NDVI數(shù)據(jù)是由最佳像元選擇合成算法處理得到的數(shù)據(jù)產(chǎn)品[7-8]。

作為實驗的驗證分析，還收集了災(zāi)中和災(zāi)后時期的GF-3和Sentinel-1觀測數(shù)據(jù)以及災(zāi)期Sentinel-2光學(xué)觀測數(shù)據(jù)和災(zāi)前Landsat影像數(shù)據(jù)。

GF-3數(shù)據(jù)來自于由中國高分應(yīng)用技術(shù)中心分發(fā)的成像方式為精細(xì)條帶1(Fine Strip I,FSI)模式，極化方式為雙極化的Level-1A級數(shù)據(jù)。

Sentinel-1和Sentinel-2數(shù)據(jù)產(chǎn)品來自于歐空局(European Space Agency，ESA)官方網(wǎng)站：(https://sci-hub.Copernicus.eu/dhus//home)的Level-1地距數(shù)據(jù)(Ground Range Detected,GRD)，成像方式為干涉寬幅(Interferometric Wide swath,IW)模式，極化方式為VV的雷達(dá)數(shù)據(jù)。Sentinel-2的數(shù)據(jù)為經(jīng)輻射校正和正射校正處理后的Level-1C級數(shù)據(jù)產(chǎn)品。

Landsat-8 OLI(Operational Land Imager)光學(xué)數(shù)據(jù)，來自于美國地質(zhì)勘查局(U.S.Geological Survey-USGS)發(fā)布的經(jīng)輻射校正和地形幾何校正處理后的L1T數(shù)據(jù)產(chǎn)品。

表1 多顆衛(wèi)星遙感觀測數(shù)據(jù)產(chǎn)品信息

二、洪災(zāi)動態(tài)信息提取方法

本研究根據(jù)上述由研究區(qū)域樣點的NDVI時序，本文提出了通過由洪災(zāi)導(dǎo)致的NDVI的變化特征來監(jiān)測和提取耕地受災(zāi)信息的方法。其主要步驟包括對最初合成產(chǎn)品中異常點的剔除以及插值處理、多時相NDVI變化特征的動態(tài)監(jiān)測、受災(zāi)程度以及災(zāi)后恢復(fù)狀況的評價。

圖2 基于多源衛(wèi)星遙感洪災(zāi)動態(tài)信息提取方法流程

(一)NDVI時序變化特征

斯里蘭卡南部主要栽培季節(jié)Maha，從9月到3月。主要灌溉的次級Yela季節(jié)從4月到9月，主要種植作物為水稻。以馬塔勒地區(qū)2017年受災(zāi)區(qū)域為例，根據(jù)當(dāng)?shù)貫?zāi)害管理中心(DMC)數(shù)據(jù)報告從中選出6個不同特征的農(nóng)田樣本區(qū)域并提取2017受災(zāi)年和2016年NDVI多時相數(shù)據(jù)。對比兩年的NDVI變化情況(圖3)。

圖3 受災(zāi)2017年(紅色)和未受災(zāi)2016年(藍(lán)色)樣點NDVI時序特征

災(zāi)前的持續(xù)強(qiáng)降水使土壤水分過度飽和以及災(zāi)后大面積農(nóng)田被水淹沒，嚴(yán)重影響農(nóng)田作物的生長，該變化會通過遙感衛(wèi)星監(jiān)測數(shù)據(jù)中的地表反射率隨時間序列的波動表現(xiàn)出來，具體表現(xiàn)為NDVI值隨洪災(zāi)持續(xù)降低。因此根據(jù)農(nóng)田的受災(zāi)變化與NDVI變化的關(guān)系特征可提取洪災(zāi)的動態(tài)信息。

1.對比未受災(zāi)年NDVI值，不同樣區(qū)受災(zāi)年NDVI值開始降低的時期基本一致

根據(jù)災(zāi)害中心(DMC)的報告得知洪災(zāi)發(fā)生的日期為5月25日。如上圖(a)(b)(c)(d)所示，受災(zāi)年NDVI值在5月末6月初下降趨勢時分明顯，表明洪災(zāi)發(fā)生后立刻對農(nóng)田作物生長產(chǎn)生了嚴(yán)重的影響。

2.不同區(qū)域受災(zāi)后NDVI值降低的持續(xù)時間不同

將兩年同時相的NDVI做差值對比發(fā)現(xiàn)受災(zāi)年同比減小約0.18，統(tǒng)計連續(xù)大于此閾值的時相數(shù)，其中(a)(b)(c)(d)分別呈現(xiàn)了8、7、4、3個時相，其中(c)(d)的NDVI值在受洪災(zāi)的影響降低，災(zāi)后又回升到了2016年的NDVI值水平。圖1(a)(b)連續(xù)有9和8個時相的NDVI差值連續(xù)大于閾值，災(zāi)后一直持續(xù)降低至年末沒有回升，該類區(qū)域表明在洪災(zāi)影響過后農(nóng)民由于嚴(yán)重的破壞而放棄后續(xù)的耕作。由此可見，連續(xù)大于閾值的時相數(shù)可以反映受災(zāi)的嚴(yán)重程度以及災(zāi)后的恢復(fù)狀況。

3.NDVI值異常降低與升高的區(qū)域

除上述受災(zāi)年NDVI值隨洪災(zāi)降低的區(qū)域之外，也存在NDVI值異常降低或升高的情況。如圖(e)從年初就一直低于為受災(zāi)年或一直高于為受災(zāi)年的變化，其因素可能為農(nóng)田變?yōu)檎晷莞蛴糜诮ㄖ玫?。如圖(f)所示異常升高又突然降低的區(qū)域可能為整個時相周期內(nèi)均為多云天氣導(dǎo)致。

(二)災(zāi)害動態(tài)信息提取

將上述耕地洪災(zāi)動態(tài)信息提取方法應(yīng)用于斯里蘭卡馬塔勒地區(qū)受到洪災(zāi)2017年和未受災(zāi)的2016年MODIS多時相NDVI數(shù)據(jù)進(jìn)行了方法驗證。首先為保證數(shù)據(jù)的準(zhǔn)確性對收集到的研究區(qū)域兩年23個時相的MODIS數(shù)據(jù)進(jìn)行異常點剔除以及線性插值等預(yù)處理；根據(jù)2017年與2016年NDVI的差值變化提取馬塔勒地區(qū)耕地洪災(zāi)動態(tài)信息。

對比兩年NDVI災(zāi)后同時相的變化差值數(shù)據(jù)表明，2017年受災(zāi)后的NDVI值比2016年同時相總體相差0.18。以此0.18為閾值，計算出NDVI降低時間Dt以及連續(xù)低于閾值的時相數(shù)Dn，根據(jù)這些特征可以評估作物受暴雨降水災(zāi)害的影響程度以及受災(zāi)后作物恢復(fù)情況。MODIS-16天NDVI產(chǎn)品數(shù)據(jù)雖然由最佳像元合成方法處理得到，很大程度上減小了云的影響[9]，但由于斯里蘭卡的熱帶季風(fēng)氣候，雨季為每年5月至8月和11月至次年2月，持續(xù)時間很長。因此會出現(xiàn)在合成周期內(nèi)均為多云天氣的狀況無法反應(yīng)樣點的真實數(shù)據(jù)，從而影響到NDVI樣點的正常變化規(guī)律。因此在提取這些特征值之前，首先對MODIS-NDVI數(shù)據(jù)存在的部分異常進(jìn)行剔除和插補。

在NDVI多時相變化中，按照農(nóng)作物的生長到成熟收割，NDVI的變化呈現(xiàn)的是逐漸增大然后減小的變化規(guī)律[10]。NDVI的變化會隨著作物在收獲季或災(zāi)害影響而減小但不會在后一時相出現(xiàn)突然增大的現(xiàn)象。因此，若樣點中連續(xù)時相之間出現(xiàn)突然增大和減小且幅度超過0.18的情況，則將其判定為受云層影響而產(chǎn)生的異常值并進(jìn)行剔除。同時結(jié)合異常點前后時相NDVI值進(jìn)行線性插值，由此減小異常值對多時相變化特征分析的影響。

將處理后得到2016年和2017年分別23個時相的NDVI數(shù)據(jù)集。通過總時相為N的受災(zāi)年(NDVIs[i])與未受災(zāi)年(NDVIc[i])相同區(qū)域與時相(i)中連續(xù)兩個時相差值ΔNDVI[i](式(1)，i=1,2,…,N)大于閾值(T)作為受災(zāi)開始時間Dt。并將其受災(zāi)后持續(xù)時相數(shù)Dn(式(2)-式(4))作為受災(zāi)持續(xù)時長。

(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)

根據(jù)上述采集到的樣點區(qū)域NDVI隨時間變化的特征結(jié)合通過閾值得出的受災(zāi)開始時間Dt以及受災(zāi)持續(xù)時間Dn，可將整個區(qū)域受災(zāi)情況劃分為4種類型。

1.災(zāi)后恢復(fù)型：受災(zāi)開始時間Dt與洪災(zāi)發(fā)生時間基本一致。同時，受災(zāi)持續(xù)時間Dn小于60天的樣點。

2.災(zāi)后棄耕型：受災(zāi)開始時間Dt與洪災(zāi)發(fā)生時間基本一致。同時，受災(zāi)持續(xù)時間Dn大于60天的樣點。

3.非災(zāi)害型：研究區(qū)域中2017年NDVI值的減小除災(zāi)害原因外，還有其他變化的因素如耕地變?yōu)榻ㄖ玫亍⑿莞氐阮愋?。根?jù)斯里蘭卡的雨季變化規(guī)律，只有在5月(即i=9)才有可能發(fā)生持續(xù)強(qiáng)降水而引發(fā)洪災(zāi)的規(guī)律，以Dt必大于9或者Dn必小于(N-9)為條件限制，可以判定如圖所示的NDVI值變化是由非災(zāi)害所引起的。

(三)結(jié)果與分析

分別為按上章所述方法逐像元處理，并按照像元分辨率得到不同特征像元的面積總數(shù)統(tǒng)計得到的受災(zāi)開始時間和持續(xù)天數(shù)的結(jié)果。圖4的開始時間是根據(jù)NDVI數(shù)據(jù)的每個時相間隔16天，將開始時相Dt進(jìn)行日期換算后的結(jié)果。圖5的持續(xù)天數(shù)也是根據(jù)每個時相間隔16天由時相個數(shù)Dn換算的天數(shù)(圖4)。

圖4 研究區(qū)受災(zāi)信息直方圖：(a)開始時期、(b)持續(xù)天數(shù)

同時，根據(jù)NDVI的變化特征以及受災(zāi)開始時間Dt、受災(zāi)持續(xù)時間Dn，根據(jù)上述不同受災(zāi)類型的劃分標(biāo)準(zhǔn)，對整個受災(zāi)區(qū)域的進(jìn)行分類提取結(jié)果(圖5)。

圖5 研究區(qū)域受災(zāi)類型提取統(tǒng)計圖

1.從受災(zāi)開始時間提取結(jié)果可以看出，大部分受災(zāi)區(qū)域的開始時間在5月24日至6月9日之間，這也與災(zāi)害中心(DMC)報道的洪災(zāi)發(fā)生時間5月25日基本一致。部分區(qū)域洪災(zāi)發(fā)生前開始降低的現(xiàn)象可能是災(zāi)害來臨之前的持續(xù)降水已經(jīng)對部分區(qū)域的作物生長產(chǎn)生了影像。

2.從受災(zāi)持續(xù)時間提取結(jié)果可以看出，大部分受災(zāi)區(qū)域影像持續(xù)時間在49至65天之間。由此可見，洪災(zāi)對大部分區(qū)域作物的影響持續(xù)了兩個月或更久，最久的區(qū)域持續(xù)了6個月。通過受災(zāi)持續(xù)時間的統(tǒng)計可以直觀的反應(yīng)不同區(qū)域受災(zāi)的嚴(yán)重程度，持續(xù)時間越長則表明受災(zāi)程度越嚴(yán)重。該數(shù)據(jù)可以為救災(zāi)計劃實施提供有力依據(jù)。

3.從受災(zāi)類型的提取結(jié)果可以看出雖有48%的受災(zāi)區(qū)域在災(zāi)后恢復(fù)到了為受災(zāi)年的水平，但仍有38%的受災(zāi)區(qū)域由于影響較為嚴(yán)重在災(zāi)后很長一段時間內(nèi)沒有恢復(fù)到年前的狀態(tài)，其中10%的非災(zāi)害類型根據(jù)聯(lián)合國糧食及農(nóng)業(yè)組織(Food and Agriculture Organization of the United Nations，F(xiàn)AO)出版的斯里蘭卡2017作物和糧食評估報告中提到由于多年的干旱導(dǎo)致斯里蘭卡糧食大幅減產(chǎn)使得部分農(nóng)田被開發(fā)為建筑用地，在高海拔地區(qū)則通過休耕改種茶樹等市場價值更高的作物。通過劃分不同的受災(zāi)類型可以直觀的反應(yīng)不同地區(qū)災(zāi)后恢復(fù)情況[10]。該數(shù)據(jù)可以為災(zāi)后損失評估提供幫助。

三、對比與驗證

SAR由于具有全天候全天時的特點，在極端天氣條件下也能夠在第一時間獲取災(zāi)區(qū)的影像信息，為救災(zāi)工作以及災(zāi)情的預(yù)警評估提供重要的依據(jù)[11]。本研究利用災(zāi)中和災(zāi)后GF-3和Sentinel-1獲取的多期數(shù)據(jù)，以更高分辨率的GF-3影像為提取對象，以同期Sentinel-2影像作為輔助進(jìn)行校正，應(yīng)用成熟的閾值分割法，提取了多期洪水淹沒范圍(圖6)。

圖6 研究區(qū)域SAR衛(wèi)星多時相水體提取結(jié)果

同時，根據(jù)不同覆蓋類型對未受災(zāi)年多時相NDVI數(shù)據(jù)進(jìn)行非監(jiān)督分類。通過地面調(diào)查和獲取到的受災(zāi)年5月28日Sentinel2光學(xué)影像(10m分辨率)的目視解譯，結(jié)合土地覆蓋分類體系不同地類NDVI多時相變化的季節(jié)性、峰值特征，分為農(nóng)田、森林、草地、灌木叢、濕地、水體、建筑用地和荒地8個類型。將上述受災(zāi)年提取的水體信息與聚類結(jié)果疊加圖(圖7)同地類被淹沒的面積隨時相的變化，由于濕地和荒地所占比例過小忽略不計。

圖7 研究區(qū)域聚類結(jié)果與各地類動態(tài)信息統(tǒng)計圖

(一)農(nóng)田淹沒面積對比驗證

由GF-3和Sentinel1衛(wèi)星數(shù)據(jù)提取到的5月30日農(nóng)田被淹沒面積為22.05km2,MODIS則取對應(yīng)于5月30日GF-3和Sentinel1時相提取的2017年與2016年的NDVI差值從該時相開始連續(xù)3個時相小大于0.18時提取的像元數(shù)。計算得出淹沒面積為25.19km2略大于GF-3衛(wèi)星提取結(jié)果，其原因可能是NDVI時序數(shù)據(jù)提取的特征像元有來自已經(jīng)退水但NDVI值無法恢復(fù)到未受災(zāi)年的區(qū)域。同時由于MODIS影像分辨率較低，無法精確提取微小水面信息，從而平滑了那些小面積的水體變化信息。

(二)農(nóng)田淹沒范圍

由GF-3和Sentinel1災(zāi)中和災(zāi)后多時相數(shù)據(jù)計算出農(nóng)田被淹沒的面積隨時相的變化可知農(nóng)田部分至8月8日水淹面積為4.23km2，相較于災(zāi)中已經(jīng)褪去82.6%。在MODIS時序中統(tǒng)計2017受災(zāi)年持續(xù)時間Dn小于65天既持續(xù)影響時間在8月8日前后的面積為19.7km2，占整個持續(xù)天數(shù)比例的78.6%，兩組數(shù)據(jù)基本一致。MODIS占比相對較少的原因可能是在水退之后由于土壤過飽和會繼續(xù)影響作物生長?？梢苑从吵龃蟛糠謪^(qū)域受災(zāi)害影響的持續(xù)時間在兩個月左右。

結(jié)束語

洪災(zāi)通過使農(nóng)田土壤水分過飽和從而影響農(nóng)作物的生長是一個漸變的過程。通過多源遙感衛(wèi)星的多時相的數(shù)據(jù)進(jìn)行動態(tài)監(jiān)測才能準(zhǔn)確評估其受災(zāi)程度和災(zāi)后恢復(fù)情況。利用獲取到的Terra/MODIS受災(zāi)年和未受災(zāi)年觀測獲取的多時相NDVI數(shù)據(jù)的變化特征，提出了通過設(shè)定閾值來提取災(zāi)害的動態(tài)信息的方法。并利用GF-3和Sentinel1衛(wèi)星對受災(zāi)區(qū)域進(jìn)行輔助，實現(xiàn)了對該區(qū)域更精確的動態(tài)監(jiān)測。結(jié)果表明，基于MODIS多時相合成的NDVI時序數(shù)據(jù)可提取不同受災(zāi)區(qū)域的開始時間、持續(xù)時間以及災(zāi)后恢復(fù)狀況等諸多時間和空間的動態(tài)信息。通過高分辨率微波遙感衛(wèi)星提取的水體信息進(jìn)行對比驗證后發(fā)現(xiàn)該方法可以快速準(zhǔn)確的獲取大范圍災(zāi)情信息。為研究暴雨、洪災(zāi)以及其他覆蓋面積較大的突發(fā)性地質(zhì)災(zāi)害提供了一個方法應(yīng)用研究案例。

通訊作者：張文君