韓曉勇，韓　玲，陳魯皖

(長(zhǎng)安大學(xué) 地質(zhì)工程與測(cè)繪學(xué)院，陜西 西安 710054)

?

基于DTW距離的時(shí)序NDVI數(shù)據(jù)植被信息提取
——以秦巴山區(qū)為例

韓曉勇，韓玲，陳魯皖

(長(zhǎng)安大學(xué) 地質(zhì)工程與測(cè)繪學(xué)院，陜西 西安 710054)

摘要：秦巴山區(qū)是我國(guó)重要的生態(tài)屏障，對(duì)該區(qū)的植被信息提取開(kāi)展研究，可為區(qū)內(nèi)生態(tài)服務(wù)功能及自然資源開(kāi)發(fā)利用提供基礎(chǔ)數(shù)據(jù)。通過(guò)加窗處理改進(jìn)DTW距離相似性算法，結(jié)合臨近度模糊分類方法對(duì)2005—2014年的MODIS NDVI時(shí)序數(shù)據(jù)進(jìn)行植被信息提取。首先利用S-G濾波對(duì)MODIS NDVI時(shí)序數(shù)據(jù)進(jìn)行重建；再利用2013年的采樣數(shù)據(jù)構(gòu)建各類植被的標(biāo)準(zhǔn)NDVI時(shí)序曲線，逐像元計(jì)算與標(biāo)準(zhǔn)NDVI時(shí)序曲線的加窗DTW距離，利用臨近度模糊分類實(shí)現(xiàn)植被信息提??；最后驗(yàn)證提取精度。結(jié)果表明，算法具有較高的運(yùn)行效率，可避免錯(cuò)誤匹配，以較高的精度(總體精度83.8%，kappa系數(shù)0.77)實(shí)現(xiàn)長(zhǎng)時(shí)間序列的植被信息提取。

關(guān)鍵詞：NDVI；時(shí)間序列；加窗DTW；臨近度；模糊分類

植被是生物圈的主要構(gòu)成部分，是陸地生態(tài)系統(tǒng)的主體，其生長(zhǎng)和分布受環(huán)境制約，可作為氣候變化的指示因子[1]。衛(wèi)星遙感從不同時(shí)間、空間尺度和不同光譜波段進(jìn)行對(duì)地觀測(cè)，獲取大量的觀測(cè)數(shù)據(jù)，在資源調(diào)查與植被變化監(jiān)測(cè)等方面得到廣泛應(yīng)用。植被指數(shù)是植被監(jiān)測(cè)的最佳指示因子之一，其中以歸一化植被指數(shù)(NDVI) 的應(yīng)用最為廣泛[2]。MODIS NDVI產(chǎn)品數(shù)據(jù)更以其高時(shí)間分辨率，較長(zhǎng)時(shí)間跨度和容易獲取的優(yōu)勢(shì)被廣泛應(yīng)用于植被信息提取及變化監(jiān)測(cè)[3-4]。不同植被類型在其生長(zhǎng)周期內(nèi)擁有不同NDVI時(shí)序特征，例如落葉林表現(xiàn)為振幅較大的單峰曲線，一年兩熟的農(nóng)田則為雙峰曲線，可利用植被自有時(shí)序特征進(jìn)行植被信息的識(shí)別[5]。Geerken利用傅里葉濾波進(jìn)行NDVI時(shí)序重建，使用相似性度量方法提取草地植被類型[6]。Wardlow通過(guò)MODIS NDVI時(shí)序決策樹(shù)來(lái)提取美國(guó)中部農(nóng)作物信息，總體精度達(dá)84%[7]。那曉東等將離散的傅里葉變換應(yīng)用于NDVI時(shí)序重建，采用傅里葉組分的相似度方法提取三江平原濕地植被信息，總體精度達(dá)79.67%，Kappa系數(shù)為0.752 5[8]。Peng利用水稻灌水期NDVI時(shí)序特征提取湖南省水稻種植面積[9]。管續(xù)棟對(duì)MODIS NDVI時(shí)序去噪后，引入動(dòng)態(tài)時(shí)間規(guī)整(DTW)距離相似性方法，提取泰國(guó)東北部地區(qū)單季稻、雙季稻面積，總體精度為83.38%[10]。

秦巴山區(qū)是我國(guó)重要的生態(tài)屏障，是南水北調(diào)中線工程的主要水源地，對(duì)該區(qū)域的植被信息實(shí)現(xiàn)高精度提取，可為其生態(tài)服務(wù)功能定位及對(duì)自然資源的合理開(kāi)發(fā)利用提供基礎(chǔ)數(shù)據(jù)。本文以2013年的野外調(diào)查數(shù)據(jù)為基礎(chǔ)，利用時(shí)序曲線相似性評(píng)價(jià)的方法提取研究區(qū)內(nèi)2005—2014年的植被特征，利用2014年的驗(yàn)證數(shù)據(jù)評(píng)價(jià)提取精度。

1研究區(qū)概況

選擇陜西省內(nèi)的秦巴山區(qū)為研究區(qū)，主要包括商洛市、漢中市和安康市，由于西安市、寶雞市南部也屬于秦巴山區(qū)，故將兩市也納入研究區(qū)(見(jiàn)圖1)，研究區(qū)總面積98 713 km2。秦嶺是黃河、長(zhǎng)江兩大水系的分水嶺和我國(guó)南北氣候分界線，北坡為暖溫帶氣候，南坡為北亞熱帶氣候，在陜西境內(nèi)連綿約500 km，海拔多在1 000 m以上，主峰太白山海拔3 767 m。巴山是嘉陵江和漢江的分水嶺，北麓為北亞熱帶氣候，南麓為中亞熱帶，海拔平均1 500~2 500 m。研究區(qū)是陜西省最大最主要的林區(qū)，屬于暖溫帶落葉闊葉林和北亞熱帶常綠落葉闊葉混交林地帶。秦嶺與巴山之間，分布著漢中盆地、西鄉(xiāng)盆地、安康盆地、商丹盆地和洛南盆地，盆地內(nèi)耕地集中，農(nóng)業(yè)生產(chǎn)水平較高，是陜南重要的農(nóng)業(yè)生產(chǎn)區(qū)和人口聚居區(qū)。此外研究區(qū)還包括部分渭河平原和渭北旱塬丘陵溝壑區(qū)。秦巴山區(qū)氣候具有明顯的過(guò)渡特征，植被類型豐富，按照國(guó)家土地利用現(xiàn)狀分類方法，結(jié)合遙感數(shù)據(jù)的分辨率特征，本文以水田、旱地(含水澆地)、有林地和灌木林4個(gè)二級(jí)類和1個(gè)草地一級(jí)類進(jìn)行特征提取。

2數(shù)據(jù)及預(yù)處理

2.1MODIS NDVI數(shù)據(jù)

從EOSDIS(http://reverb.echo.nasa.gov/)下載研究區(qū)2005—2014年的MODIS植被指數(shù)產(chǎn)品(MOD13A2)，空間分辨率1 km，16 d一期數(shù)據(jù)，10年共230期數(shù)據(jù)。使用MRT工具將每期HDF文件中NDVI和NDVI 質(zhì)量控制文件進(jìn)行拼接、裁剪，統(tǒng)一轉(zhuǎn)換到WGS84 UTM N49投影帶，并利用質(zhì)量控制文件剔除質(zhì)量差的像元。

2.2野外采樣與驗(yàn)證數(shù)據(jù)

2013年4月至7月及9月對(duì)研究區(qū)內(nèi)植被情況進(jìn)行野外調(diào)查，共調(diào)查樣點(diǎn)173個(gè)，為彌補(bǔ)野外調(diào)查樣點(diǎn)的不足，本文利用2013年、2014年World view的高清影像，結(jié)合野外實(shí)測(cè)數(shù)據(jù)的判讀特征確定130個(gè)精度驗(yàn)證點(diǎn)，兩類樣點(diǎn)的空間分布見(jiàn)圖1。

圖1　研究區(qū)地理位置及樣點(diǎn)數(shù)據(jù)分布

3研究方法

3.1S-G濾波算法原理

理論上植被的NDVI時(shí)序曲線應(yīng)是連續(xù)平滑的，但由于傳感器自身性能、太陽(yáng)光照角度、觀測(cè)視角及云、大氣氣溶膠等觀測(cè)條件和隨機(jī)干擾因素的影響，導(dǎo)致NDVI時(shí)序曲線呈鋸齒狀的不規(guī)則波動(dòng)，反映季節(jié)變化趨勢(shì)不明顯，限制NDVI時(shí)序在植被覆蓋變化分析和信息提取等領(lǐng)域的應(yīng)用，為此發(fā)展一系列NDVI時(shí)序重建的方法[11]。選用較常用的Savitzky-Golay(S-G)濾波方法對(duì)MODIS NDVI產(chǎn)品進(jìn)行數(shù)據(jù)重建，S-G濾波是一種移動(dòng)窗口的加權(quán)平均算法，其設(shè)計(jì)思想是利用n階多項(xiàng)式來(lái)擬合滑動(dòng)窗口內(nèi)時(shí)序數(shù)據(jù)，多項(xiàng)式系數(shù)使用最小二乘法得出，NDVI時(shí)序數(shù)據(jù)集的S-G濾波過(guò)程可描述為

(1)

式中： Y*為S-G時(shí)序重建數(shù)據(jù);Yj+i代表原始時(shí)序數(shù)據(jù);Ci為濾波系數(shù);N為歸一化因子;m為滑動(dòng)窗口寬度。

通過(guò)不同參數(shù)設(shè)置對(duì)比，選用二階多項(xiàng)式，滑動(dòng)窗口寬度(m=2)為5，歸一化因子N為35，這既保證擬合效果，又避免過(guò)度擬合。由于滑動(dòng)窗口寬度為5，濾波處理后重建的時(shí)序的長(zhǎng)度減少為19，即損失最初兩期(1月份)和最后兩期(12月份)的數(shù)據(jù)，這兩個(gè)月份對(duì)于植被特征標(biāo)志意義不大，故本文只用濾波后每年19期重建數(shù)據(jù)進(jìn)行相似性計(jì)算。

3.2加窗DTW距離

利用2013年調(diào)查數(shù)據(jù)的NDVI均值曲線作為各植被類型標(biāo)準(zhǔn)生長(zhǎng)曲線，再逐年計(jì)算各像元時(shí)序曲線與標(biāo)準(zhǔn)生長(zhǎng)曲線的DTW距離。年際間植被生長(zhǎng)的水熱條件、光照條件不盡相同，從而導(dǎo)致植被每年的NDVI時(shí)序曲線發(fā)生不同程度的扭曲，因此DTW距離更適合該類時(shí)序曲線相似性評(píng)價(jià)。

(2)

(3)

邊界條件：w1=dmn，wK=d11，該路徑的起止元素為距離矩陣的斜對(duì)角線的兩端元素；連續(xù)性：若ws=dr,c，ws-1=dr′c′，0≤r-r′≤1 且0≤c-c′≤1，保證路徑中相鄰元素的連續(xù)性；有界性：max(m,n) ≤K≤m+n-1，路徑所經(jīng)過(guò)的矩陣元素個(gè)數(shù)K存在上限和下限。

3.3臨近度模糊分類

從模糊分類的角度，遙感數(shù)據(jù)的混合像元即不完全屬于或完全不屬于某個(gè)類別，使用貼近度來(lái)描述模糊集與標(biāo)準(zhǔn)模糊集靠近程度，將混合像元貼近度最高的類別作為該像元的歸屬類別。模糊集定義為待識(shí)別的像元與各植被類型參考時(shí)序曲線的DTW距離，標(biāo)準(zhǔn)模糊集是各類別參考時(shí)序曲線間的DTW距離，式(4)~式(6)為貼近度計(jì)算過(guò)程。

(4)

(5)

(6)

4結(jié)果及精度評(píng)價(jià)4.1標(biāo)準(zhǔn)時(shí)序曲線

使用迭代處理消除樣點(diǎn)數(shù)據(jù)受混合像元及植被類型變化的影響。首先提取S-G 濾波后的NDVI樣點(diǎn)數(shù)據(jù)，逐期計(jì)算各類均值作為初始參考曲線，再計(jì)算各樣點(diǎn)與初始參考曲線的DTW距離，將DTW距離大于2倍標(biāo)準(zhǔn)差的樣點(diǎn)進(jìn)行剔除，將剩余樣本點(diǎn)計(jì)算新的參考曲線，重復(fù)迭代過(guò)程，迭代計(jì)算的控制條件是前后兩次參考曲線的DTW距離變化小于0.000 1，最終確定各類標(biāo)準(zhǔn)時(shí)序曲線(見(jiàn)圖2)。圖2中水田和旱地NDVI值的年內(nèi)變化均為雙峰結(jié)構(gòu)，這與一年兩熟的生產(chǎn)方式一致；水田的NDVI上升時(shí)間最早，結(jié)束較晚，生長(zhǎng)期比旱地區(qū)長(zhǎng)。草地、灌木林與有林地均為單峰曲線，65—129日三者的NDVI增長(zhǎng)較快，此時(shí)為展葉期到快速生長(zhǎng)階段；129—273日NDVI處于平穩(wěn)高值區(qū)，植被穩(wěn)定生長(zhǎng)；273日以后NDVI快速衰落，植被進(jìn)入休眠階段。

圖2　各地類標(biāo)準(zhǔn)時(shí)序曲線

4.2加窗DTW距離算法效率

加窗DTW距離算法的設(shè)計(jì)目的是減小計(jì)算量，本文使用IDL語(yǔ)言將傳統(tǒng)DTW算法和加窗

DTW算法進(jìn)行對(duì)比測(cè)試，在CPU主頻3.3 G HZ、內(nèi)存4 G的硬件配置下，計(jì)算研究區(qū)383行518列每年19期數(shù)據(jù)的平均耗時(shí)分別為：傳統(tǒng)DTW算法42.815 73 s，加窗DTW算法25.744 16 s，加窗DTW算法可提高39.87%的計(jì)算效率。另外，加窗算法能有效避免錯(cuò)誤匹配，本文對(duì)水田和旱地樣點(diǎn)數(shù)據(jù)進(jìn)行錯(cuò)誤匹配測(cè)試，加窗算法可避免18.6%的錯(cuò)誤匹配。

4.3模糊分類

表1為各類參考時(shí)序曲線之間的DTW距離，表中的每列(或每行)是對(duì)應(yīng)地類的標(biāo)準(zhǔn)模糊集?？梢?jiàn)水田、旱地與草地、有林地、灌木林的DTW距離較大，說(shuō)明水田、旱地與另外3類之間較易區(qū)分。灌木林與草地和有林地的DTW距離較小，說(shuō)明三者之間易發(fā)生錯(cuò)分，尤其是灌木林與有林地的DTW距離為0.013 8，二者之間錯(cuò)分的可能性最大。

表1　標(biāo)準(zhǔn)模糊集

按照貼近度最大原則對(duì)研究區(qū)2005—2014年的植被類型進(jìn)行特征提取，結(jié)果見(jiàn)圖3。圖中各地類信息的空間分布合理，水田主要分布在漢中、西鄉(xiāng)和安康盆地，旱地分布在商丹和洛南盆地以及渭河平原區(qū)，有林地主要分布在秦嶺與巴山山區(qū)，灌木林和草地主要分布在耕地與有林地的過(guò)渡地帶。表2為各植被類型面積統(tǒng)計(jì)，2005—2014年有林地面積緩慢增加了約2 500 km2，旱地面積減少較多，尤其在2011年陜西實(shí)施的“避災(zāi)扶貧移民搬遷工程”后旱地面積顯著減少，減少的主要是山區(qū)的坡耕地，其它3類的面積變化較小。

表2　各植被類型面積統(tǒng)計(jì)　km2

4.4分類精度

使用130個(gè)驗(yàn)證點(diǎn)對(duì)提取結(jié)果進(jìn)行精度分析，得到分類混淆矩陣見(jiàn)表3。本文算法總體分類精度為83.8%，Kappa系數(shù)0.77，說(shuō)明利用DTW距離與模糊分類相結(jié)合的方法可以較高精度提取各類植被。水田與旱地的制圖精度較低，二者的漏分誤差較高，圖2中旱地與水田的分布特征也顯示該特點(diǎn)。草地與灌木林的用戶精度較低，說(shuō)明兩類的錯(cuò)

表3　分類精度

分誤差較大，圖3中可見(jiàn)灌木林與草地主要分布在有林地與耕地的過(guò)渡地帶，這類地物在遙感影像中主要以混合像元形式存在，從而導(dǎo)致它們的用戶精度較低。

利用2006—2014年的《陜西省統(tǒng)計(jì)年鑒》中各市(區(qū))農(nóng)作物播種面積對(duì)旱地和水田分類結(jié)果進(jìn)行檢驗(yàn)，旱地和水田面積的平均相對(duì)誤差分別為15.3%和14.5%，最大相對(duì)誤差分別為23.8%和 19.7%，分類精度整體較高。

5結(jié)論

本文使用2005—2014年的NDVI數(shù)據(jù)進(jìn)行植被信息提取，時(shí)間跨度較大，基于DTW距離的相似性分析有效解決年度間氣像因子差異引起的曲線偏移，并通過(guò)加窗處理有效提高DTW算法效率，減少錯(cuò)誤匹配，結(jié)合臨近度模糊分類能夠避免直接選取閾值造成的因植被物候期差異造成的植被誤分類，且植被信息提取結(jié)果總體精度較高。

然而在長(zhǎng)時(shí)序提取過(guò)程中仍然存在一些問(wèn)題有待進(jìn)一步研究。首先提取結(jié)果受到MODIS NDVI產(chǎn)品數(shù)據(jù)質(zhì)量的制約。例如，研究區(qū)內(nèi)2009年樣點(diǎn)的NDVI數(shù)據(jù)S-G濾波后，DTW距離全部超過(guò)2倍標(biāo)準(zhǔn)差，最終2009年的信息提取結(jié)果嚴(yán)重失真，這應(yīng)該是與2009年存在11期大面積異常數(shù)據(jù)有關(guān)，這意味著若年度數(shù)據(jù)存在近一半的噪聲數(shù)據(jù)時(shí)要設(shè)計(jì)其他方法才能保證特征提取的質(zhì)量；其次是混合像元影響。混合像元的類別只代表像元內(nèi)分布最廣的植被類型，例如在西安東北部土地利用方式多為旱地與果園混雜分布，而果園的NDVI特征與有林地相近，從而產(chǎn)生在某些年度有林地-旱地之間反復(fù)的情況，從而降低信息提取精度，可通過(guò)結(jié)合作物物候和中高分辨率遙感影像進(jìn)行混合像元分解，以提高提取精度實(shí)現(xiàn)植被類型進(jìn)一步細(xì)分。

參考文獻(xiàn)：

[1]任園園,張哲,侯欽磊,等.中國(guó)大陸植被覆蓋對(duì)降水與溫度變化的響應(yīng)[J].水土保持學(xué)報(bào)，2013,33(2)：78-82.

[2]陳燕麗，羅永明，莫偉華，等.MODIS NDVI 與 MODIS EVI 對(duì)氣候因子響應(yīng)差異[J].自然資源學(xué)報(bào)，2014,29(10):1802-1809.

[3]陳燕麗，龍步菊，潘學(xué)標(biāo),等.基于MODIS NDVI和氣候信息的草原植被變化監(jiān)測(cè)[J].應(yīng)用氣象學(xué)報(bào)，2010,21(2):229-236.

[4]楊斌，王金生.基于GIS的丘陵區(qū)耕地景觀格局時(shí)空演變特征分析[J].測(cè)繪工程，2014,23(9):1-4.

[5]李文葉,姜魯光,李鵬.2001—2010年鄱陽(yáng)湖圩區(qū)水稻多熟種植時(shí)空格局變化[J].資源科學(xué)，2014，36(4)：809-816.

[6]GEERKEN R,ZAITCHIK B,EVANS J P.Classifying rangeland vegetation type and coverage from NDVI time series using Fourier Filtered Cycle Similarity[J].International Journal of Remote Sensing,2005,26(24):5535-5554.

[7]WARDLOW B D,EGBERT S L,KASTENS J H.Analysis of time-series MODIS 250m vegetation index data for crop classification in the US Central Great Plains[J].Remote Sensing of Environment,2007,108(3):290-310.

[8]那曉東,張樹(shù)清,李曉峰,等.MODIS NDVI時(shí)間序列在三江平原濕地植被信息提取中的應(yīng)用[J].濕地科學(xué),2007,5(3):227-236.

[9]PENG D，HUETE A R，HUANG J,et al.Detection and estimation of mixed paddy rice cropping patterns with MODIS data[J].International Journal of Applied Earth Observation and Geoinformation，2011，13(1)：13-23.

[10] 管續(xù)棟,黃翀,劉高煥,等.基于DTW距離的時(shí)序相似性方法提取水稻遙感信息：以泰國(guó)為例[J].資源科學(xué),2014,36(2):267-272.

[11] 李杭燕,頡耀文,馬明國(guó).時(shí)序NDVI數(shù)據(jù)集重建方法評(píng)價(jià)與實(shí)例研究[J].遙感技術(shù)與應(yīng)用，2009，24(5):596-602.

[責(zé)任編輯：張德福]

Extraction of vegetation information using adding windows DTW distance with NDVI time series data

HAN Xiaoyong,HAN Ling，CHEN Luwan

(School of Geology Engineering and Geomatics，Chang’an University，Xi’an 710054，China)

Abstract：Qinling-Bashan Mountain area is an ecological barrier.The vegetation information extraction research can provide base data for ecological service function and natural resources exploitation.DTW distance similarity measure is improved by adding windows，and vegetation information is extracted by proximity-fuzzy classification from MODIS NDVI in the years of 2005 to 2014.S-G filter method is first applied to weakening the noise and reconstructing the NDVI time series.Each type of vegetation standard NDVI time series curves is generated from the 2013 sample data.For each to-be-classified pixel，a quantitative similarity between its time-curve and standard time series curves is calculated using adding windows DTW，then vegetation information is extracted in research area.Finally，validate data is used to verify the extraction accuracy.The result shows that the algorithm has high efficient and can avoid mismatch，and can realize long time serial vegetation information extraction with high accuracy (overall accuracy:83.8%，Kappa coefficient:0.77).

Key words：NDVI；time series；adding windows DTW；proximity；fuzzy classification

中圖分類號(hào)：TP79

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1006-7949(2016)03-0011-06

作者簡(jiǎn)介：韓曉勇(1981-)，男，博士研究生.

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目(41201099)

收稿日期：2015-06-11