基于紋理特征與MODIS-NDVI時間序列的耕地面積提取研究*

于利峰，烏蘭吐雅※，烏云德吉，許洪滔，包珺瑋，任婷婷

(1.內(nèi)蒙古自治區(qū)農(nóng)牧業(yè)科學院農(nóng)牧業(yè)經(jīng)濟與信息研究所，呼和浩特 010031； 2.內(nèi)蒙古自治區(qū)農(nóng)業(yè)遙感工程技術研究中心，呼和浩特 010031)

0 引言

隨著國家快速發(fā)展，耕地保護的狀態(tài)和形勢不容樂觀，耕地面臨大量被占用的壓力。中國是一個擁有十幾億人口的農(nóng)業(yè)大國，耕地是糧食安全的重要前提與保障，也是實施國家“藏糧于地”重要科技戰(zhàn)略的基礎注http：//www.gov.cn/zhengce/content/2017-04/10/content_5184613.htm。因此，及時了解耕地變化情況，掌握耕地利用信息，對于制定相關土地政策有重要意義。

歸一化植被指數(shù)(Normalized Difference Vegetation Index，NDVI)是國際上主要用于識別植被的重要參數(shù)[1]，它與植被蓋度有明顯的相關性[2]，在研究區(qū)域和全球的土地覆蓋變化、植被生長過程等方面有廣泛應用。同一區(qū)域、時間序列相同的植被，其NDVI具有相似的變化曲線。利用這一特點可以通過時間序列方式來研究耕地和其他土地利用類型的NDVI變化方式，并以此提取耕地面積。

1999年，美國航空航天局的中分辨率成像光譜儀(EOS-MODIS)成功發(fā)射，該數(shù)據(jù)的光譜、時間分辨率高，其NDVI值對植被的變化更敏感[3]，是目前研究中最常用的數(shù)據(jù)。在利用MODIS-NDVI時序數(shù)列的分類研究中，邊金虎[4]等采用若爾蓋高原濕地MODIS-NDVI時間序列產(chǎn)品，分析了濕地分布情況； 司亞輝[5]等，采用最大值合成法，得到MODIS-NDVI的每旬合成數(shù)據(jù)，并結合地面數(shù)據(jù)推斷錫林郭勒盟的草地類型在空間布局情況； 宮攀[6]利用MODIS旬最大值合成數(shù)據(jù)，結合物候參數(shù)，提高了土地利用分類精度。但是，地物類型的物候參數(shù)有時候很類似，MODIS-NDVI時序數(shù)列有可能表現(xiàn)出誤導性的地物變化特征。基于此文章將紋理特征加入到MODIS-NDVI時間序列數(shù)據(jù)集中參與分類，并對其在耕地面積提取中表現(xiàn)進行比較、分析。該文將對進一步使用MODIS-NDVI數(shù)據(jù)開展土地利用分類提供了理論依據(jù)。

圖1 研究區(qū)地理位置

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

研究區(qū)地處內(nèi)蒙古呼倫貝爾市，位于119°28′E～120°34′E、49°06′N～49°28′N，該區(qū)域雨熱同期，蒸發(fā)量大，夏季多云(圖1)。研究區(qū)深受海洋季風和潮濕氣流的影響，由東向西逐漸減弱，依次形成濕潤區(qū)、半濕潤區(qū)和半干旱區(qū)[7]，相應的地帶性植被從東向西為森林、草甸草原和草原。該區(qū)域地塊面積大，分布著機械化程度高、勞動生產(chǎn)率高以及商品率高的國有農(nóng)場及私人農(nóng)場，耕地連片集中，適合使用遙感監(jiān)測[8]。

1.2 MODIS-NDVI數(shù)據(jù)處理

該文使用的數(shù)據(jù)為2016年MODIS13Q1，空間分辨率為250m，時間分辨率為16d，選擇h25v3、h25v4、h26v3、h26v4的共92景影像，在NASA網(wǎng)站獲取。利用MRT(MODIS Reproject Tools)軟件對數(shù)據(jù)進行鑲嵌和格式轉換，并將投影由原來的正弦投影轉換為Albers_WGS_1984。最后在ENVI中對像元異常值進行剔除，對于小于-1的像元值賦值為-1，大于1的像元值賦值為1，并對去除異常值的數(shù)據(jù)波段合成，形成一組具有23個波段多波段時序數(shù)據(jù)。

1.3 MODIS-NDVI數(shù)據(jù)時間序列重構

研究區(qū)多云氣候使得16d合成的MODIS-NDVI數(shù)據(jù)仍然存在云影噪聲，需要對其濾波降噪。Savitzky-Golay平滑濾波器(以下簡稱S-G濾波)通常用于平滑頻率跨度較大的噪聲數(shù)據(jù)[9]。最初由Savitzky和Golay于1964年提出，發(fā)表于Analytical Chemistry 雜志[10]。在數(shù)據(jù)流的光滑降噪等方面被廣泛使用，是一種在固定時間范圍內(nèi)根據(jù)多項式擬合的濾波手段。該濾波方法最大的特征是可以在降低噪聲的同時可以保證信號的形狀、寬度不變[11]。當信號的典型峰值窄時，該濾波器工作特別好，通常會保留曲線的高度和寬度，算法如式(1)：

(1)

式中xk.smooth、xk+i分別是濾波擬合結果與原始數(shù)組，k是原數(shù)組的角標。m是某點所在右側或左側包含的數(shù)據(jù)點個數(shù)即濾波的窗口大小，m的值越大，會產(chǎn)生了更平滑的結果，代價是平坦化尖峰，值通常取3～7。ci是濾波器的系數(shù)，該系數(shù)由m的高階多項式的階數(shù)決定，這里可以將階數(shù)看作為D，它的典型值為2～4，階數(shù)較低會產(chǎn)生更平滑的結果，但可能引入偏差； 階數(shù)較高將降低濾波器偏差，但可能會過度擬合數(shù)據(jù)并產(chǎn)生較嘈雜的結果，并且階數(shù)要小于2m+1[12]。該文選擇的數(shù)據(jù)是16d合成數(shù)據(jù)，S-G濾波能有效去除云對NDVI數(shù)據(jù)的影響，使結果更傾向于真實數(shù)值(圖2)。

鑒于不同的濾波參數(shù)所獲得的濾波效果不同，該文將參數(shù)分別設置為m=3，D=2、m=5，D=3進行結果比較。(圖3)

圖2 S-G濾波去除云效果

圖3 研究區(qū)各地表類型S-G平滑前后NDVI時間序列對比

圖3中耕地、林地、草地、居民地、未利用土地和水域之間NDVI時間變化曲線區(qū)別比較明顯，林地的NDVI相對于其他地物一直處于較高水平，在它的變化曲線中161～225d之間時間序列曲線趨之平緩，NDVI值在0.8～0.9之間，符合北方5月下旬到8月下旬林木生長的旺盛情況； 草地的NDVI時間變化曲線比較平緩，最大值在0.5左右，在97～161d、225～289d草地的NDVI值均大于耕地，這段區(qū)間分別是耕地的作物出苗期和成熟期，草地的長勢略好于耕地。

耕地在129d開始呈幾何式增長，到177d達到峰值，此時正值作物出苗到長成時期，植被指數(shù)不斷增高； 從209～273d是作物植被指數(shù)第一個下降時期，此時研究區(qū)作物陸續(xù)進入成熟期，作物生長情況步入衰退期； 273d后是作物植被指數(shù)第二個下降時期，這個階段研究區(qū)作物基本收割完畢，研究區(qū)天氣在開始變得寒冷，作物收割后的裸地或者干枯的秸稈是耕地的主要背景。

在非植被地物類型中，3種地物類型NDVI值均處于較低水平，水體的主要特征是NDVI一直為0以下； 居民地與未利用地時間序列曲線變化趨勢類似，但未利用地的在65～305d區(qū)間要略高與居民地，可以用作分類的特征比較。

由對比圖3可以看出，隨著擬合窗口大小(m)與多項式擬合階數(shù)(D)的增大，時間序列曲線更平滑，擬合窗口大小等于5，擬合階數(shù)等于3時，出現(xiàn)了明顯的數(shù)據(jù)失真與過度擬合情況。因此該文在測試對比后，選擇m=3，D=2作為Savitzky-Golay的濾波參數(shù)，借助于IDL編程，對MODIS-NDVI數(shù)據(jù)時間序列進行重構。

1.4 紋理特征提取

紋理由灰度像素在空間位置上反復顯現(xiàn)而形成的[13]，因此圖像空間中相隔一段位置的兩個像素之間存在相對的灰度聯(lián)系，即圖像中灰度的空間相關特性[14]。灰度共生矩陣則是用來描述紋理的常用手段，該矩陣是20世紀70年代初由R.Haralick等人總結的，是目前認為的在提取紋理特征時一種有效的統(tǒng)計方法[15]。研究區(qū)位于內(nèi)蒙古呼倫貝爾市，該地區(qū)耕地集中連片，大面積的耕地表面存在紋理，且特征明顯。通過反復測試，選擇3×3窗口，對其中的6個統(tǒng)計量進行提取，式(2)～(7)：

(2)

Variance=∑i∑j(i-u)2p(i，j)

(3)

(4)

Entropy=-∑i∑jp(i，j)log(p(i，j))

(5)

SecondMoment=∑i∑j{p(i，j)}2

(6)

(7)

式中，n為灰度值的階數(shù)，p(i，j)是n×n的歸一化共生矩陣，u為p(i，j)的均值，σxσy為p(i，j)的相關系數(shù)。如圖4是以第161d的NDVI數(shù)據(jù)為例，計算6個統(tǒng)計量的結果，可以看到耕地所在區(qū)域灰度特征明顯，紋理特征集中。

圖4 第161d的NDVI數(shù)據(jù)灰度共生矩陣統(tǒng)計量計算結果

圖5 主成分分析獲得紋理特征參數(shù)的主成分分布

圖6 實測驗證數(shù)據(jù)分布

數(shù)據(jù)冗余往往會導致處理效率變低，在對NDVI時間序列的23個波段計算灰度共生矩陣時會產(chǎn)生大量數(shù)據(jù)，因此采用PCA對時間序列23個波段紋理特征共138統(tǒng)計量進行分析，提取主要成分。通過主成分分析，發(fā)現(xiàn)前7個主要成分對影像的紋理解釋的累積貢獻率為95.97%，因此取前7個主成分參與之后的分類。

2 研究方法

2.1 分類體系的選擇

該文主要研究內(nèi)容是耕地信息的提取，在分類體系選擇時，主要依據(jù)中科院中國1： 10萬比例尺土地利用現(xiàn)狀遙感監(jiān)測數(shù)據(jù)庫的一級分類標準，包括耕地、林地、草地、水域、居民地和未利用土地6個類型[16]。對于大面積宏觀監(jiān)測來說，若要通過普查地物來選擇訓練樣本十分困難，因此在建立感興趣區(qū)時，參照部分2016年呼倫貝爾市部分無云Landsat8數(shù)據(jù)、谷歌地球數(shù)據(jù)和實測點選取分類樣本。(圖6)

通過上述方法，共選取訓練樣本3 506個像元(表1)，選擇其中的60%參與分類，剩余40%進行結果驗證。

表1 樣點地物類型及數(shù)量統(tǒng)計

樣點地物類型像元個數(shù)占比(%)耕地53315.202林地81723.303草地82223.445水域62417.798居民地3349.526未利用地37610.724

2.2 分類器的選擇

決策樹分類是利用遙感數(shù)據(jù)及其他輔助數(shù)據(jù)，利用知識經(jīng)驗總結、統(tǒng)計和數(shù)學歸納方法等，獲得提取規(guī)則建立決策樹并進行數(shù)據(jù)分類[17]。Breiman等于1984年提出的分類與回歸樹(Classification and Regression Tree，CART)是一種通過對測試數(shù)據(jù)和目標數(shù)據(jù)構成的決策數(shù)據(jù)集構建的決策樹結構[18]。該算法在研究區(qū)有足夠的樣本數(shù)量的情況下，采用GiNi系數(shù)作為最佳測試變量和分割閾值的準則[19]。基尼系數(shù)定義如式(8)(9)：

(8)

(9)

式(9)中，P(j/h)是從訓練數(shù)據(jù)集里隨機選取的一個樣本；nj(h)是訓練數(shù)據(jù)里測試變量值為h時屬于第j類的樣本數(shù)；n(h)是訓練數(shù)據(jù)里該檢測變量值為h的個數(shù)；j為類別數(shù)。

圖7 技術流程

2.3 研究路線

在構建MODIS-NDVI時間序列數(shù)據(jù)后，使用Savitzky-Golay濾波器對時間序列數(shù)據(jù)進行平滑降噪處理，并計算紋理特征，通過CART決策樹算法對不同數(shù)據(jù)提取耕地、林地、草地、水域、居民地和未利用土地6個類型地物面積，具體流程如圖7。

3 結果與分析

3.1 可分離性評價

Jeffries-Matusita參數(shù)常用于表示2個類別間的差異程度[20]，這兩個參數(shù)在提取時應大于1.85； 如果小于1.8，需要重新選擇樣本； 如果參數(shù)值比1小，則應該將兩類合并為一類。通過地面實測點數(shù)據(jù)，選擇一定數(shù)量的訓練樣本，對原始的MODIS-NDVI時間序列、重構的MODIS-NDVI時間序列以及添加紋理信息濾波后的MODIS-NDVI時間序列3種數(shù)據(jù)計算各種土地類型之間的可分離度，均達到了1.95以上屬于合格樣本，可以用于分類。

3.2 分類結果及對比分析

利用S-G濾波方法對MODIS-NDVI時間序列數(shù)據(jù)進行重構，地物類型之間通過時間序列曲線的差異用于分類提取：林地是各類地物中生長季比較長的類型，研究區(qū)的林地NDVI在129～161d之間增長變化較快，且能夠在161～225d內(nèi)保持較高的NDVI 值，這與其他類別差別較大； 草地的NDVI時間序列比耕地的NDVI時間序列寬，二者的NDVI最大值也都在0.5左右，草地低于耕地的NDVI最大值。基于此特征，利用CART方法進行土地利用分類提取耕地信息。如圖8，依次為原始MODIS-NDVI時間序列、重構MODIS-NDVI時間序列以及紋理+重構MODIS-NDVI時間序列的分類結果。

從分類結果的簡單分析可以發(fā)現(xiàn)，分類結果總體符合預期要求，由于MODIS-NDVI數(shù)據(jù)的空間分辨率為250m，很多細節(jié)沒有被體現(xiàn)，比如河流、村級居民點沒有被識別出來。從宏觀結果來看，3種分類結果所體現(xiàn)出的土地利用相似，耕地基本分布在大興安嶺沿麓地區(qū)，包括莫力達瓦自治旗、阿榮旗、扎蘭屯市、鄂溫克自治旗及額爾古納市，與真實情況相同。從細節(jié)上看，3種分類結果仍有不同之處，如圖8。

通過與高分辨率影像對照發(fā)現(xiàn)，MODIS數(shù)據(jù)由于分辨率較低，分類結果會產(chǎn)生許多破碎像元，可以通過聚類濾波將小圖斑剔除。S-G濾波后MODIS-NDVI時間序列更加精確地區(qū)分地物，能夠將耕地特征準確識別出來，在圖中使用原始MODIS-NDVI時間序列分類結果中許多耕地信息沒有被表現(xiàn)出來(圖9a、圖9b)。而加入紋理信息后，使得地物的表面結構組織排列屬性展現(xiàn)得更加豐富，圖9呼倫湖在添加紋理特征后湖岸界限分明，表明紋理特征對于大面積連篇地物的區(qū)分作用更為明顯(圖9c、圖9d)。上述結果雖然仍存在錯分小圖斑，但已基本滿足宏觀監(jiān)測需求。

3.3 分類結果的精度驗證

選擇之前未參與分類的實測點與高分辨率數(shù)據(jù)訓練的感興趣區(qū)合計1 564個像元，對分類結果進行評價，具體Kappa系數(shù)和分類的精度見表2～4。

表2 原始數(shù)據(jù)的分類精度與Kappa系數(shù)

地物類型林地耕地草地未利用地水體居民地制圖精度用戶精度林地17613500099.4490.72耕地01941301877.689.81草地1434559273678.0471.77未利用地00107137274157.8143.91水域0002142079.7898.61居民地003615438.8584.38Kappa系數(shù)=0.662 總體精度=74.04%

表3 S-G濾波的分類精度與Kappa系數(shù)

地物類型林地耕地草地未利用地水體居民地制圖精度用戶精度林地17762001010086.76耕地0204006181.696.68草地04054711452793.8374.62未利用地0016116315748.9552.73水域0003135075.8497.83居民地000405438.8593.1Kappa系數(shù)=0.719 總體精度=78.84%

圖8 數(shù)據(jù)分類

圖9 分類結果細節(jié)對比

表4 S-G濾波結合紋理特征的分類精度與Kappa系數(shù)

地物類型林地耕地草地未利用地水體居民地制圖精度用戶精度林地17762001010086.76耕地0206107182.495.81草地0385446542693.3180.35未利用地0018167156370.4663.5水域0001151084.8399.34居民地000404935.2592.45Kappa系數(shù)=0.789 總體精度=83.72%

表5 耕地面積提取準確度對比

數(shù)據(jù)提取出耕地面積(萬hm2)準確率(%)原始NDVI時間序列249.98245.61重構NDVI時間序列194.54779.85紋理+重構NDVI時間序列184.05586.33

從表2～4顯示結果可以看到，分類精度在通過對NDVI時間序列數(shù)據(jù)濾波平滑和添加紋理特征后總體分類精度提高了約10%。其中林地的分類精度最高，主要原因是其地物特征明顯，而且遙感數(shù)據(jù)上的NDVI時間變化特征最為典型。此外未利用地分類精度最差，其與草地、城鎮(zhèn)地物類型的混淆錯誤是導致總體分類精度下降的原因，從時間序列曲線中不難發(fā)現(xiàn)草地與未利用地在春、秋、冬3個時間段趨勢與值幾乎相同； 未利用地與居民地的NDVI變化趨勢近似。未利用地與居民地、草地之間區(qū)分度不夠明顯是導致錯分、誤分的主要原因。

3種分類結果的分類精度從高到低依次為S-G濾波結合紋理信息的分類結果、單獨S-G濾波的分類結果、原始NDVI時間序列結果。其中，S-G濾波結合紋理信息的分類結果總體精度達到了83.72%，Kappa系數(shù)為0.789，滿足了大尺度土地利用制圖要求。

該文主要目的是提取耕地信息。為此，在查閱2016年內(nèi)蒙古自治區(qū)統(tǒng)計年鑒可知，研究區(qū)耕地面積為161.925萬hm2，通過比較耕地面積的準確度，來評價耕地信息提取的優(yōu)劣。具體的提取出的耕地面積與準確度見表5。

通過以上分析得出：NDVI作為一種特征參數(shù)，其所構成的時間序列可以滿足大尺度的耕地信息提取研究，而重構后的NDVI時間序列能更加準確描述地表的各類地物變化情況，利于特征提取，特別是與紋理信息結合之后分類精度有明顯提高，說明紋理特征有效增強像素之間相對的灰度聯(lián)系，描述了圖像中灰度的空間相關特性，是遙感分類中的關鍵信息，該文還為今后對此類信息的進一步探索提供了參考。

4 結論

該文以MODIS-NDVI時間序列為數(shù)據(jù)源，采用Savitzky-Golay濾波器對NDVI時間序列進行重構，并在此基礎上結合主成分分析后的紋理特征信息用于土地利用分類和耕地信息提取。通過結果對比，得到以下結論：

(1)紋理特征作為該文的主要輔助信息，彌補了單一NDVI時間序列上存在的不足，擴大了NDVI數(shù)據(jù)像元間的像元差異性，添加紋理特征后的分類結果精度顯著提升，總體分類精度為83.72%，其中耕地面積提取的準確度達到了86.33%，對于耕地資源大面積調(diào)查具有積極意義。

(2)MODIS-NDVI時間序列各地物的曲線表現(xiàn)出的變化特征是分類的重要指標，其中林地、水域特征較為明顯，居民地與未利用地間差異較小，因而林地與水域類型在MODIS-NDVI時間序列地表分類中有明顯的優(yōu)勢。

(3)該文通過對MODIS-NDVI時間序列濾波處理提高了分類精度，經(jīng)過Savitzky-Golay濾波后能更真實反映地物的NDVI時間變化特征，而紋理信息在一定程度上提高了分類的精度和準確度，添加紋理信息的數(shù)據(jù)能夠綜合提升約10%的土地利用分類精度，增加40.72%的耕地提取準確度。

(4)在文中的精度驗證時使用了真實地表感興趣區(qū)對分類結果進行了評價，雖然感興趣不能完全反映整體分類精度結果，還有可能導致結果虛高，但筆者認為均勻分布的真實地表感興趣區(qū)仍可以代表一定程度的評價結果。

下一步該文將在數(shù)據(jù)選取、處理與研究方法上加以補充，進一步提高MODIS數(shù)據(jù)在大尺度農(nóng)業(yè)資源監(jiān)測中的優(yōu)勢。