利用NDVI時間序列識別汶川地震滑坡的分布L1楊文濤;汪明;史培軍【摘要】In2008,theWenchuanearthquakecausedmassivelandslidesinmountainousregions.Earthquakeinducedlandslidescandistributewidelyandcanhavespecialspatialdistributionpatterns.Themainlimitsofusinghighresolutionsatelliteimagestoidentifylandslidesareeasyaffectionbyclouds,highcostsandlongerreturnperiod.Toovercomethesefaults,itisnecessarytoexploreaneffectivemethodtoidentifylandslidesbyusingthefreecostandhightemporalresolutionMODISdata.ByusingMODISNDVIproductsof16dayinterval,yearlyNDVItimeseriescanbeformedtoidentifylandslidesafterSavitzky-Go-layfiltering.ThespatialdistributionoflandslidesinducedbythiseventcanbefurtheridentifiedbycomparingNDVIseriesbetween2007and2008aswellasbyconsideringtheoccurringtimeofthegreatearthquake.Inthisresearch,Pingwucountyisselectedasanexperimentalareatotesttheefficiencyoftheproposedmethod,andlandslidesinterpretedfromETMimageswith30mresolutionareusedtotestitsprecision.BasedonaboveanalysisandbyaddingslopeinformationfromDEM,thismethodisthenappliedtothewholeWenchuanaffectedregionwhichhas13county-leveldistrictstoproducelandslidedistributionmapsofthisregion.%利用MODISNDVI產品，生成了全年16天間隔的NDVI時間序列,再經過Savitzky-Golay濾波后消除了云的影響，可以用來識別滑坡.假設2007年、2008年研究區土地利用狀況基本不變，對比2007年和2008年兩年的NDVI序列,并結合汶川地震的發生時間，可以識別出研究區發生滑坡的區域分布.結合平武縣30mETM目視解譯滑坡結果，對自動識別方法進行了驗證結果表明:MODIS識別出的滑坡與目視解譯結果匹配較好.在此基礎上，將本文方法應用到受汶川地震影響最嚴重的13個縣區級行政區，并綜合研究區坡度信息，識別出了汶川重災區滑坡分布概圖.【期刊名稱】《遙感信息》【年(卷)，期】2012(027)006【總頁數】5頁(P45-48,56)【關鍵詞】汶川地震;NIVI時間序列;Savitzky-Golay濾波;滑坡識別;滑坡分布【作者】楊文濤;汪明;史培軍【作者單位】民政部-教育部減災與應急管理研究院，北京100875;北京師范大學地表過程與資源生態國家重點實驗室，北京100875;民政部-教育部減災與應急管理研究院,北京100875;北京師范大學地表過程與資源生態國家重點實驗室，北京100875;北京師范大學環境演變與自然災害教育部重點實驗室，北京100875;民政部-教育部減災與應急管理研究院，北京100875【正文語種】中文【中圖分類】TP791引言滑坡災害，通常是由外界因素觸發不穩定坡體引起的，地震和強降雨是引發滑坡的兩種最常見的觸發因子。山地地區尤其容易發生滑坡災害，而這些地區又多被濃密的植被覆蓋。在滑坡識別工作中，遙感技術是最常用的手段。但是，由于這些覆蓋有茂密植被的山地地形往往起伏很大、山體陰影以及云的影響非常顯著，使得遙感對滑坡的識別不能很好地發揮作用，而后續的觸發因子(余震或降水)會繼續作用，從而形成多次觸發因子下的滑坡分布。建立單次滑坡觸發因子下的完整滑坡分布對地質災害風險評估至關重要。傳統的滑坡解譯技術主要利用航空立體像對識別滑坡，但是這種方法對影像的時相要求較高，需要盡快在單次觸發滑坡的因素發生之后立即獲得影像；與此同時，航空遙感成本很高，對天氣條件要求也高，這些都是利用航空遙感監測滑坡分布的障礙［1］。隨著QuickBird,IKONOS,SPOT等衛星的出現，高分、超高分辨率遙感衛星影像成為滑坡目視識別的重要信息源［2~3］。此外，利用滑坡前后高分辨率DEM的差別，滑坡空間分布以及滑坡體的體積也可以得到識別［4］。但是，高分、超高分衛星影像具有時間分辨率低、重訪周期長、易受云影響等缺點。利用衛星影像識別滑坡，主要是對比滑坡前后同一地區的影像，總體上分為兩大類方法：基于像元尺度的滑坡體識別方法和基于影像形態、分割、模式識別的方法。—般而言，非基于像元的方法優于基于像元的方法［5］。隨著LiDAR(LightDetectionAndRanging)技術的引入，單個滑坡體可以得到精確的量測［6~7］。與此同時，應用PS(PermanentScatters)技術，干涉雷達(InSAR)也可以用于監測地表毫米級的微小形變［8］。美國NASA的MODIS傳感器安裝在Aqua和Terra衛星平臺上，空間分辨率有250m,500m,1km等幾種，粗空間分辨率帶來的是重訪周期短、低成本等優勢。文章采用MODIS13Q1數據產品，該數據空間分辨率為250m，時間分辨率為16天。2008年5月12日，四川省汶川縣發生里氏8.0級大地震，震中位于31.0°N,103.4°E［9］。這次大地震引發了超過56000個滑坡，覆蓋了41750km2的面積［10］，這些滑坡多發在㈣度以上列度區［11］；在隨后的幾個月里又發生了無數次余震，這些余震引發了一系列后續的地震型滑坡。同時，震后4個月里災區不斷發生強降雨事件，又進一步引發了地震-降雨型滑坡和降雨型滑坡。滑坡分布面積較廣，且當時的天氣狀況較差，給滑坡的遙感調查、解譯帶來了一定的難度。文章旨在利用低成本的MODIS粗分辨率衛星影像，研究滑坡事件的空間分布。2研究區概況與數據研究區包括四川、甘肅地震重災區的13個市縣區，地理坐標位于：103°50'E~104°58'E、31°59'N~33°02'N之間，屬青藏高原到四川盆地的過渡地帶，境內多高山，海拔從600多米到5000米以上，變化很大。研究區地勢西北高東南低，地形崎嶇，山區深切溝谷發育，地質構造復雜，地質災害多樣且頻發。本文使用的MODIS數據是2007年和2008年的MODIS13Q1數據，該數據是以16天為周期，空間分辨率為250m的NDVIL3產品。驗證數據是30m分辨率的Landsat7ETM影像兩景，覆蓋平武縣東部地區，影像獲取時間是2008年4月30日和2008年5月16日。將2007年和2008年的各23張NDVI產品分別生成23波段的〃多波段”數據，即影像中每個像元有23個波段，代表從當年1月到12月16天間隔的NDVI時間序列。植被的生長狀況具有年際循環周期的物候特征［12］:春季不斷加速生長，生物量開始增加；夏季加速生長，葉綠素含量、生物量達到峰值；秋季開始逐漸下降；冬季降到最低點。由于MODIS的NDVI產品容易受到云的影響，并且常常在反映植被狀況的時候偏低。云對NDVI的影響具有突變性，即云對NDVI值有降低的影響。但滑坡對植被擾動造成的NDVI的影響則具有迅速下降和緩慢恢復的表現。針對云對NDVI的影響，陳晉等人［13］在2004年給出了Savitzky-Golay濾波算法，通過該算法，可以有效地去除云對NDVI時間序列的影響，使原始NDVI像元值更加接近于真實值。Savitzky-Golay濾波方法的原理是采用最小二乘卷積法平滑NDVI時間序列。這種卷積窗口可以看作是一個權重窗口，作用是消除低值（被認為是噪音弓I起的），保留高值。這種濾波適用于連續或等間隔的信號序列處理。平滑的濾波公式為：其中，Y是原始的NDVI值，是濾波之后的NDVI值，Ci是窗口中第i個NDVI值的系數或權重，N是窗口的大小，j是NDVI時間序列中第j個需要處理的NDVI值。平滑濾波效果是由窗口的大小和濾波多項式的級數決定的：窗口越小，濾波多項式級數越低，平滑效果越不明顯，細節保存越多；反之，細節保存較少。文章選擇濾波多項式級數為6，半窗口大小為4。圖1自動識別流程3滑坡識別該方法建立在兩個假設的基礎上：研究區主要被植被覆蓋；植被的變化由滑坡造成，即選取的研究區大面積被植被覆蓋且在2007年到2008年之間土地利用類型沒有發生變化。在逐像元計算時，單像元在2007年與2008年的時間序列相似；如果時間序列有變化，即在2008年5月發生明顯的下降，在2008年后期出現恢復，則判斷該時間序列所屬像元是汶川地震引發的滑坡，識別該像元為滑坡點位。結合上述像元時間序列判別標準，首先對平武縣境內滑坡自動識別，流程見圖1,共識別出了全縣5.9km2的滑坡面積，滑坡發生區主要集中在平武縣東南沿涪江的谷地及平武縣東部汶川斷裂帶附近，其他地區滑坡分布較少，如圖2所示。驗證數據為兩期Landsat7ETM數據，驗證區位于平武縣東部滑坡分布較為密集的地方，總面積約13.8km2（圖3）。圖2MODIS識別的滑坡分布圖3滑坡識別驗證（a）2008年4月30日ETM影像；（b）文獻［14］中研究區地質災害三維視圖；（c）文獻［14］中采用的2008年5月16日2.5m分辨率SPOT5解譯的滑坡；（d）目視解譯的滑坡體分布與MODIS滑坡識別結果的對比通過對比汶川地震前后驗證區的影像，同時參考了《5.12汶川地震典型地質災害影像研究》［14］中平武縣地區地質災害三維影像圖，目視解譯出了滑坡區域范圍（其中，圖3最大的滑坡的是平武縣馬鞍山滑坡［15］）。在此基礎上，把MODISNDVI識別滑坡與目視解譯的滑坡區域疊加（圖3）。對ETM目視解譯結果與MODIS自動識別結果定量對比（表1）。結果表明：MODIS方法識別出的滑坡面積占ETM識別滑坡總面積的90.79%，二者完全匹配的比例為55.02%,有44.98%的目視解譯滑坡沒有被識別出來，35.77%的面積發生誤判。從圖3以及整體解譯結果來看，能被MODISNDVI自動正確識別的滑坡體一般面積較大，基本為能完整覆蓋一個MODIS像素的滑坡體。其中，大于4個MODIS像元的滑坡其識別精度最高。而研究區較多的不能被識別的滑坡多是由于面積較小而不完全覆蓋整個或大部分MODIS單像素的小型滑坡造成的。這些誤差往往是由于MODIS像素覆蓋下的地表景觀，非滑坡性因素的信息占主導地位所致，這也和采用MODIS數據較低的空間分辨率有著直接的關系。表1ETM目視解譯與MODIS滑坡識別比較在保證精度的基礎上，將研究區擴大至汶川地震重災區，包括：甘肅省隴南市武都區、文縣，四川省青川縣、平武縣、江油市、北川縣、安縣、茂縣、綿竹縣、什那縣、彭州市、都江堰市、汶川縣、理縣等14個縣市區，包含了全部的IX度及以上地震烈度區，大部㈣度區，以及部分VD度區和小部VI度區。圖4汶川地震重災區滑坡識別分布圖通過提出的NDVI自動識別技術，得到了汶川重災區滑坡解譯分布。滑坡的發生與坡度有很大關系，在原算法的基礎上，考慮了地形坡度對滑坡發生的影響，即對坡度小于5°的地形不予考慮滑坡的發生。根據這一條件，最后得到了汶川地震重災區滑坡識別分布圖（圖4）。從結果圖來看，汶川重災區滑坡主要沿東北-西南向狹長帶狀分布，這與龍門山中央斷裂帶走向基本一致。4結束語由滑坡造成的植被的擾動可以反映在MODISNDVI像元值的變化上，進而反映在整年的NDVI時間序列變化上。對比正常年與地震年份同一像元NDVI時間序列的變化，可以識別出植被擾動發生的區域，進而識別出滑坡發生的區域。利用MODIS的16天NDVI時間序列產品識別之前，需要經過Savitzky-Golay濾波，以消除時間序列中的噪音和由于云的影響造成的NDVI的降低。將通過NDVI時間序列提取出的滑坡分布，和ETM目視解譯的結果進行了驗證，結果表明：MODISNDVI時間序列識別出的滑坡與較高分辨率衛星解譯的結果有—定的一致性，可以利用較粗分辨率的MODIS產品識別滑坡發生的空間分布，然后有針對性地選取滑坡多發區購買高分、超高分影像識別滑坡，從而快速、準確地建立完整的滑坡分布。美國MODIS衛星的重訪時間達每天2次，可以充分利用這種高時間分辨率的特性，建立間隔更短的NDVI時間序列分辨率來識別滑坡。結合正常植被的物候特征，利用MODISNDVI時間序列識別不同時間階段滑坡的分布，這在滑坡體動態變化研究中具有重要的意義。同時，可以利用單一事件（如震后的某次強降雨）訓練特定滑坡區域的NDVI序列特征，來提高預測大區域滑坡空間分布的準確性。盡管可以利用MODIS時間序列，快速、連續地監測大面積分布的滑坡災害，但是災區土地利用變化對受災區滑坡監測帶來的影響也不容忽視。因此，在滑坡災害的自動檢測研究中，利用災害發生前后最近時間的衛星影像是十分必要的，這樣可以避免災后土地利用變化對滑坡識別造成的誤差，在這一點上，MODIS的較高時間分辨率有一定的優勢。參考文獻[1]WESTENCJ,CASTELLANOSE,KURIAKOSESL.Spatialdataforlandslidesusceptibility,hazard,andvulnerabilityassessment:Anoverview[J].EngineeringGeology,2008(102):112-131.[2]MANTOVANIF,SOETERSR,WESTENV.RemotesensingtechniquesforlandslidestudiesandhazardzonationinEurope[J].Geomorphology,1996,15(3-4):213-225.[3]METTERNICHTG,HURNIL,GOGUR.Remotesensingoflandslides:Ananalysisofthepotentialcontributiontogeo-spatialsystemsforhazardassessmentinmountainousenvironments[J].RemoteSensingofEnvironment,2005,98(2-3):284-303.[4]DEWITTEO,DemoulinA.MorphometryandkinematicsoflandslidesinferredfrompreciseDTMsinWestBelgium[J].NaturalHazardsandEarthSystemSciences,2005(5):259-265.[5]MARTINJB,FRANKLINS.DetectingtranslationallandslidescarsusingsegmentationofLandsatETM+andDEMdatainthenorthernCascadeMountains,BritishColumbia[J].RemoteSensing,2003,29(40):510-517.[6]JONESLD,WHITWORTHK.LandslideLaserScanning:Anewlookatanoldproblem[J].QuarterlyJournalofEngineeringGeologyandHydrogeology,2003(36):155-157.[7]JONESL.MonitoringlandslidesinhazardousterrainusingterrestrialLiDAR:AnexamplefromMontserrat[J].QuarterlyJournalofEngineeringGeologyandHydrogeology,2006,39(4):371-373.[8]LINH,ZHAOQ,JIANGLM,etal.Astudyofgrounddeformation