基于MODIS－NDVI時序數據的黃土丘陵區土地覆蓋分類研究

陸廣勇，楊勤科2，，王海江

（1.中國科學院 水利部 水土保持研究所，陜西 楊凌712100；2.西北大學，西安710069）

土地覆蓋是全球變化研究的重要內容，其研究為區域和全球尺度的氣候和生態過程模型提供了重要的參數。區域尺度的土地覆蓋信息主要基于多光譜遙感影像Landsat MSS／TM／ETM＋數據提取，但所需數據量龐大，而且關鍵生長時期的影像會受到云的影響而無法獲取，數據更新受限［1］。

土地覆蓋變化具有季節性，這種季節變化可以利用植被指數的變化進行定量分析研究［2］。衛星遙感影像能夠監測植被的物候變化［3－5］，來自遙感影像的區域植被指數NDVI值隨時間周期性升高和降低是植被生長周期的典型體現［6］。反映地表覆蓋物季節性變化特征的多時相NDVI可以作為土地覆蓋的分類特征［5］。

植被指數時間序列的數據源主要包括NOAAAVHRR、SPOT－VGT和MODIS等植被指數數據產品。前兩者的合成周期為10d，空間分辨率約為1 km，多用于全球或大區域的研究［3，5，7－8］。MODIS植被指數產品合成周期為16d，但其空間分辨率較高，有500m 和250m 兩種分辨率產品［9－11］，因此應用研究也越來越多［1，12－16］。

基于植被指數時間序列影像的土地覆蓋分類方法可分成非監督分類、監督分類和分層決策樹分類3種。盛永偉［3］、李道峰［7］、宮攀［1］等基于氣象衛星的NDVI時間序列影像、運用非監督分類法進行了植被覆蓋類型分類；劉勇洪等基于MODIS地表反射率、植被指數和地表溫度數據產品，選擇6種分類特征，通過監督分類，將山東省土地覆蓋分成了8類［13］；張霞等綜合運用ISODATA非監督分類法、決策樹法，對華北平原進行了土地覆蓋分類［14］。

但是很多研究應用中，土地覆蓋的季節性變化規律沒有被作為主要的分類特征，而且非監督分類法與決策樹法的應用中都是在動態聚類的基礎上進行的，物理意義不夠明確。目前在地形復雜區（如黃土丘陵溝壑區）基于MODIS植被指數時間序列影像的土地覆蓋分類研究較少。本實驗綜合運用MODIS地表反射率與NDVI植被指數兩種數據產品、DEM，結合當地物候特征，分析了土地覆蓋植被指數內變化特征，運用植被指數、濕度指數和紋理特征進行分層決策樹分類，利用ISODATA法和當地土地覆蓋分布特征調整分層決策樹結果，為區域環境動態監測、土壤侵蝕評價提供更高時間分辨率的基礎數據。

1 研究區域概況

延河流域地處陜北黃土高原中部，地理位置為36°23′－37°17′N 和 108°45′－110°28′E，面 積 達7 684km2。該流域屬暖溫帶大陸性半干旱季風氣候，年平均氣溫9℃，年均降雨量為500mm，其中7－9月降雨量占全年降雨量的60%以上。流域植被區劃上屬于森林草原地帶。由于土壤侵蝕嚴重，土壤類型以黃土母質上發育的黃綿土為主，土壤質地均一，土質疏松，抗侵蝕能力差。流域內溝壑縱橫，地形復雜，植被覆蓋率低，屬典型黃土丘陵溝壑區［17－18］。

2 研究方法

2.1 數據源

（1）本文采用NASA USGS提供的植被指數數據產品MOD13Q1的歸一化植被指數MODIS－NDVI，空間分辨率為250m，合成周期為16d。從2006年1月1日至2006年12月30日，共23個時相影像。

（2）NASA USGS提供的地表反射率產品數據MOD09GA，空間分辨為500m，影像獲取時間為2006年9月9日，共7個波段，但第5波段存在壞線，故本實驗使用了MODIS的1－4，6－7共6個波段的影像。

（3）2006年9月10日成像的Landsat TM 5圖像解譯生成的2006年土地利用圖（簡稱LU2006）和矢量化的2000年延河流域土地利用圖（簡稱LU2000）。

圖1 延河流域地理位置及地勢圖

（4）DEM：利用1∶5萬地形圖和ANUDEM軟件生成的10m分辨率DEM，用以生成坡度圖。

2.2 數據預處理

（1）基于時間序列諧波分析法（HANTS）對MODIS－NDVI進行濾波。MODIS－NDVI影像局部容易受到云等因素的影響，需要進行平滑處理［15］。本文選擇了HANTS對23景MODIS－NDVI進行平滑處理。該方法的核心主要是通過傅里葉變換和最小二乘法，利用諧波重構時間序列數據［19］。通過HANTS濾波后，植被指數時間序列更加光滑（見圖2），更有利于有效揭示物候周期性變化規律，便于地表覆蓋信息提取［19－20］。

圖2 平滑處理前后典型土地覆蓋類型NDVI時間序列曲線

（2）延河流域植被年內變化可劃分為休眠階段、恢復階段、生長季節、退化階段4個時段［18］。基于IDL語言編程，計算平滑后全年和生長期（GD，growth duration）MODIS－NDVI時間序列的最大值（max、GDmax）、最小值（min、GDmin）、差值（最大值－最小值，D、GDd）。（3）對平滑后的全年 MODIS－NDVI影像進行主成分分析（PCA）。基于反映年內植被覆蓋差異較大耕地信息PCA第二組分，在3×3窗口下，進行紋理特征局部平滑（Homogeneity，Homo），以便提取耕地信息。

（4）利用ENVI系統的GLT（Geometry Lookup Table）對 MOD09GA和 MODIS－NDVI進行幾何校正［21］。將濾波后的 MODIS－NDVI及其生成數據、MOD09GA進行Albers等面積投影。并以DEM為基準，首先將TM 5影像和2006年土地利用圖等進行空間配準，然后基于配準后的TM 5影像對MODIS－NDVI、MOD09GA進行配準。空間配準的誤差都控制在一個像元以內。

（5）基于 MOD09GA提取了影像的濕度指數（NDMI，Normalized Difference Water Index）［13］，其表達式如式（1）。

式中：ρ4，ρ6——經過標準化處理的 MODIS可見光綠色波段（第4波段）和短紅外波段（第6波段）。波段ρ6處于水的吸收帶，對濕度、含水量信息尤為敏感，而ρ4對水體的反射率較高，因此利用NDMI提取水體的外邊界。

2.3 影像分類方法

不同的土地覆蓋類型具有不同的植被時間序列曲線。根據當地的年內植被變化特征，基于植被指數響應曲線特征，進行影像分類。本文綜合運用MODIS數據的時間、光譜、空間信息，選取了植被指數NDVI、濕度指數、紋理特征作為分層決策樹的分類特征，基于閾值進行分層決策樹分類，然后運用ISODATA法分類別進行動態聚類，調整決策樹分類結果，最后利用DEM調整耕地類型，最終形成土地覆蓋類型圖。

（1）分析NDVI時間序列曲線，參考LU2000、LU2006和植被年內變化特點［18］，選擇各個指標閾值，進行分層決策樹分類，并生成各土地覆蓋類型掩膜。如圖3所示，耕地的NDVI曲線在生長期內（第113～289d）上升和下降梯度最大，波峰明顯，差值最大；典型林地的NDVI值全年都較大，在生長期內最小值略小于0.4；草地也存在較明顯的波峰，但其變化梯度小于耕地，且峰值一般小于林地和耕地；而居民地和水體NDVI年內變化較小，因為居民地內往往與一些植被（如城市綠地）混合，所以居民地的NDVI值略大于水體。根據各種土地覆蓋的時間響應曲線特征，選擇合適的閾值，進行決策樹分類。

（2）基于決策樹分類類型掩膜，利用迭代自組織數據分析法（ISODATA法），對 MODIS－NDVI圖像進行動態聚類（非監督分類）［14］，參考土地利用圖和MOD09GA多光譜影像，調整決策樹分類結果。MODIS－NDVI主要反應植被信息，因此，基于MODIS－NDVI時間序列影像對林地、高覆蓋草地進行非監督分類。在植被蓋度較低區域，NDVI指標反應植被蓋度不夠靈敏，單純依靠NDVI時間序列難以正確識別土地覆蓋類型［6］，因此基于地表反射率影像MOD09GA，對低覆蓋草地、水體與居民地進行非監督分類。

圖3 決策分類樹示意圖

（3）退耕還林（草）實施后，延河流域內大于25°的坡耕地基本都進行了退耕，因此基于DEM形成的坡度圖，將以上遙感分類結果中大于25°的耕地調整為草地。考慮到所用數據源分辨率比較低，因而本次研究將要劃分的土地利用類型為耕地、林地、草地、居民地和水體5類。

3 結果分析

3.1 空間格局與數量統計分析

基于上述分類方法得到延河流域2006年土地覆蓋分類結果，結果見圖4。植被分布由東南到西北逐漸減少，與流域年降雨量空間分布一致，林地主要分布在流域的南部，與李正國的研究結果一致［18］，占流域總面積的15.9%；草地遍布整個流域，是主要土地覆蓋類，占流域總面67.3%，其中約15%的草地是根據坡度土從遙感影像分類結果中的耕地調整而來；耕地主要沿溝道分布，在流域上游也呈現集中分布，為流域總面積的16.5%；居民地與水體兩者所占比重不足1%，主要原因是兩種流域面積比重較少，分布零散，在250m的影像上，像元混合現象嚴重［14］。

圖4 2006年土地分類結果

3.2 精度分析

選擇2006年基于TM的土地利用／覆蓋分類圖LU 2006進行比較，同時參考LU 2000，進行分類精度評價。精度分析在ENVI 4.7軟件下完成，隨機選擇耕地、草地和林地各覆蓋類型的0.1%作為驗證像元；由于居民地、水體的像元較少，故每一類各選擇10個像元作為驗證，5類地表覆蓋類型共152個像元，占總流域面積近0.12%［22］。混淆矩陣（見表1）表明總體分類精度為75.00%，Kappa系數為0.620 9，除了居民地，其它4類的用戶精度到大于70%，這就說明分類效果良好。用戶精度是分類結果可靠性指 標［6］，是 最 常 用 的 精 度 指 標［1，13－14，23］。 表 1表明，實際地表的林地和耕地被錯分為草地的情況較為嚴重，分類結果中草地類型被高估，耕地和林地被低估。這主要原因在于延河流域地表破碎，林地、耕地分布不集中，在影像分辨率較粗的情況下，像元內信息趨于平滑。同時，居民地與水體兩者地表反射率都較低，反射特征相近，所以兩者之間也存在混分現象。居民地（包括交通用地）和水體的斑塊一般較小，在粗分辨的遙感影像上很難識別，因此分類結果中兩者容易被低估。

表1 誤差矩陣

4 結論與討論

（1）以季節性變化特征作為分類基礎，根據植被指數時序曲線、影像紋理和濕度指數等分類特征，進行分層決策樹分類，并結合其它輔助數據（如地表反射率、DEM等），按類別進行非監督分類，最后調整土地覆蓋分類結果獲得土地覆蓋分類圖，具有一定的物理意義。分類結果與2006年土地利用總體一致性為75.00%，Kappa系數為0.620 9，圖像分類結果良好，能為區域環境監測、土壤侵蝕評價提供可靠的年際土地覆蓋基礎數據。

（2）根據退耕還林后大于25°陡坡上耕地稀少的特點，利用DEM生成的坡度圖，調整遙感分類結果中的耕地類型，能在一定程度上減小土地覆蓋類型分布破碎對分類結果帶來的影響，改善分類效果。

（3）MODIS－NDVI雖然在空間分辨率上較NOAA－AVHRR和SPOT－VGT有了提高，但是在地形復雜、斑塊破碎的黃土丘陵溝壑區，在250m的影像上像元混合現象嚴重，小斑塊地類較難識別，導致居民地、水體在分類結果中面積比重小于實際水平。

（4）分類結果中相當一部分耕地和林地被錯分為草地，主要原因是當地土地覆蓋類型分布破碎，影像分辨率較粗，像元混合現象嚴重。下一步工作可以結合當地植被的物候特征、分布特點，通過NDVI時間序列曲線的變化特征進行更加細致地定量分析，參考混合像元分解方法，對區域土地覆蓋類型進行分類，提高分類精度。

［1］ 宮攀，陳仲新，唐華俊，等.基于 MODIS溫度／植被指數的東北地區土地覆蓋分類［J］.農業工程學報，2006，22（9）：94－99.

［2］ 布和敖斯爾，劉紀遠，吳祖南.基于季相及經度特征的中國土地覆蓋變化遙感研究［J］.地理學報，1998，53（12）：52－60.

［3］ 盛永偉，陳維英，肖乾廣，等.利用氣象衛星植被指數進行我國植被的宏觀分類［J］.科學通報，1995，40（1）：68－71.

［4］ Jonsson P，Eklundh L.Seasonality extraction by function fitting to time－series of satellite sensor data［J］.Ieee Transactions on Geoscience and Remote Sensing，2002，40（8）：1824－1832.

［5］ 汪權方，李家永.基于時序NDVI數據的中國紅壤丘陵區土地覆被分類研究［J］.農業工程學報，2005，21（2）：72－77.

［6］ 趙英時.遙感應用分析原理與方法［M］.北京：科學出版社，2003.

［7］ 李道峰，田英，郝芳華.基于NDVI數據的黃河流域地表植被覆蓋綜合分類研究［J］.水土保持研究，2003，10（4）：88－91.

［8］ 盧玲，李新，董慶罕，等.SPOT4－VEGETATION中國西北地區土地覆蓋制圖與驗證［J］.遙感學報，2003，7（3）：214－221.

［9］ Vegetation Indices 16－Day L3Global 500m［EB／OL］.https：／／lpdaac.usgs.gov／lpdaac／products／modis＿products＿table／vegetation＿indices／16＿day＿13＿global＿500m／mod13a1.

［10］ Vegetation Indices 16－Day L3Global 250m［EB／OL］.https：／／lpdaac.usgs.gov／lpdaac／products／modis ＿products＿table／vegetation＿indices／16＿day＿l3＿global＿250m／myd13q1.

［11］ 王正興，劉闖，Huetealfredo.植被指數研究進展：從AVHRR－NDVI到 MODIS－EVI［J］.生態學報，2003，23（5）：979－987.

［12］ Guerschman J P，Hill M J，Renzullo L J，et al.Estimating fractional cover of photosynthetic vegetation，nonphotosynthetic vegetation and bare soil in the Australian tropical savanna region upscaling the EO－1Hyperion and MODIS sensors［J］.Remote Sensing of Environment，2009，113（5）：928－945.

［13］ 劉勇洪，牛錚.基于 MODIS遙感數據的宏觀土地覆蓋特征分類方法與精度分析研究［J］.遙感技術與應用，2004，19（4）：217－224.

［14］ 張霞，孫睿，張兵，等.基于 MODIS植被指數時間譜的華北平原土地覆蓋分類［J］.農業工程學報，2006，22（12）：128－132.

［15］ 左麗君，張增祥，董婷婷，等.MODIS／NDVI和 MODIS／EVI在耕地信息提取中的應用及對比分析［J］.農業工程學報，2008，24（3）：167－172.

［16］ 王紅說，黃敬峰.基于 MODIS NDVI時間序列的植被覆蓋變化特征研究［J］.浙江大學學報：農業與生命科學版，2009，35（1）：105－110.

［17］ 謝紅霞.延河流域土壤侵蝕時空變化及水土保持環境效應評價研究［D］.西安：陜西師范大學，2008.

［18］ 李正國，王仰麟，吳健生.等.不同土地利用方式對黃土高原植被覆蓋季節變化的影響：以陜北延河流域為例［J］.第四紀研究，2005（06）：762－769.

［19］ Roerink G M，Menent W.Verhoef.Reconstructing cloudfree NDVI composites using Fourier analysis of time series［J］.International Journal of Remote Sensing，2000，21（9）：1911－1917.

［20］ 王丹，姜小光，唐伶俐，等.利用時間序列傅立葉分析重構無云 NDVI圖像［J］.國土資源遙感，2005（2）：29－32.

［21］ 張樹譽，王釗，李星敏.提高衛星遙感資料利用率的方法［J］.測繪學院學報，2004，21（1）：30－33.

［22］ 劉詠梅，楊勤科，溫仲明.地形復雜地區遙感圖像分類方法應用研究：以黃土丘陵溝壑地區坡耕地遙感調查為例［J］.水土保持通報，2003，23（4）：30－32.

［23］ Kumar U N，Kerle T V.Ramachandra.Constrained Linear Spectral Unmixing Technique for Regional Land Cover Mapping Using MODIS Data［J］.Innovations and Advanced Techniques in Systems，Computing Sciences and Software Engineering，2008：416－423，DOI：10.1007／978－1－4020－8735－6＿78.