基于MODIS數(shù)據(jù)的漳滏河灌區(qū)作物種植結(jié)構(gòu)提取

高葉鵬， 劉海新， 王冬利， 文韶鑫， 錢以臨

(河北工程大學礦業(yè)與測繪工程學院，河北邯鄲 056038)

作物種植結(jié)構(gòu)作為農(nóng)業(yè)生產(chǎn)活動結(jié)果的直接體現(xiàn)，可以反映農(nóng)作物的種植類型、種植布局、種植面積以及種植模式的綜合農(nóng)業(yè)信息[1]?？焖?、準確地獲取農(nóng)作物種植結(jié)構(gòu)的時空動態(tài)變化信息，是研究作物長勢、優(yōu)化種植結(jié)構(gòu)的重要依據(jù)，也是指導農(nóng)業(yè)生產(chǎn)、估測農(nóng)作物產(chǎn)量、實施農(nóng)作物精細化管理的重要支撐[2]。因此，開展作物種植結(jié)構(gòu)提取研究具有非常重要的現(xiàn)實意義。目前，農(nóng)作物種植結(jié)構(gòu)提取的研究主要考慮遙感數(shù)據(jù)源的選擇、特征變量的選取以及分類方法的選用3個方面[3]。在數(shù)據(jù)源的選擇方面，目前常被學者用于研究的是高時間、低空間分辨率的MODIS和中高空間分辨率、時間分辨率較低的Sentinel-2、Landsat、高分系列等數(shù)據(jù)[4-6]。在特征變量選取方面，常被國內(nèi)外學者用來研究農(nóng)作物分類的特征變量是紋理特征、光譜特征、植被指數(shù)等，其中植被指數(shù)作為表征地表植被生長狀況的重要參數(shù)應用最廣泛[7-13]。在分類方法選用方面，國內(nèi)外學者較常用且分類精度較高的是隨機森林、決策樹、面向?qū)ο蟮确诸惙椒?，這些方法都較成熟，已經(jīng)得到廣大學者的驗證與應用[14-17]。但是根據(jù)已有研究結(jié)果[18-20]，這些分類方法都存在以下問題：一是需要人工采集一定數(shù)量的訓練樣本；二是需要過多的人為干擾；三是自動化程度較低。因此，本研究基于遙感技術(shù)在種植結(jié)構(gòu)領域的應用，選用MODIS數(shù)據(jù)作為數(shù)據(jù)源，以增強植被指數(shù)(EVI)作為主要特征變量，分析基于決策樹和廣義DEM模型分類方法的河北省邯鄲市漳滏河灌區(qū)主要作物種植結(jié)構(gòu)提取，對比評價2種分類方法的提取精度和分類結(jié)果，分析討論廣義DEM在提取作物方面的適用性以及提取優(yōu)勢。同時，獲得區(qū)域內(nèi)農(nóng)作物的種植結(jié)構(gòu)信息，以期為漳滏河灌區(qū)作物的精細化管理和可持續(xù)發(fā)展提供科學依據(jù)。

1 研究區(qū)概況與數(shù)據(jù)處理

1.1 研究區(qū)概況

1.1.1 地理位置 漳滏河灌區(qū)位于河北省邯鄲市中部，介于36°15′～36°55′N、113°50′～115°20′E之間，范圍涉及磁縣、永年區(qū)、曲周縣、雞澤縣、臨漳縣、成安縣、魏縣、廣平縣、肥鄉(xiāng)區(qū)、館陶縣、大名縣、邯鄲市區(qū)(邯山區(qū)、叢臺區(qū))、峰峰礦區(qū)，是全國特大型灌區(qū)之一，灌區(qū)內(nèi)總耕地面積達到26多萬hm2。灌區(qū)內(nèi)主要取水模式為水庫引水，大多來源于滏陽河上的東武仕水庫和漳河上的岳城水庫，另灌區(qū)內(nèi)建設有總干渠、干渠、分干渠等多個灌溉渠道，可實現(xiàn)2個水庫之間的聯(lián)合調(diào)度、水源共享(圖1)。

1.1.2 農(nóng)作物種植情況 漳滏河灌區(qū)屬于溫帶半干旱、半濕潤季風氣候，氣候條件適宜小麥、玉米、辣椒等農(nóng)作物的種植以及蘋果、梨、桃等林果的生長，且該地區(qū)冬小麥和夏玉米大多采取連種的方式，其主要農(nóng)作物物候期見表1。

表1 主要農(nóng)作物物候信息

1.2 數(shù)據(jù)來源及處理

1.2.1 遙感數(shù)據(jù) 本研究采用從美國國家航空航天局(National Aeronautics and Space Administration，NASA)官網(wǎng)(https：//ladsweb.modaps.eosdis.nasa.gov/search/)獲取的2018年空間分辨率為250 m，16 d合成的MOD13Q1影像，共計23個時相作為時間序列數(shù)據(jù)集，并利用MRT(MODIS Reprojection Tool)對其進行拼接、重投影、波段提取，通過裁剪以獲取研究區(qū)增強植被指數(shù)(EVI)、紅光波段反射率、近紅外波段反射率等多個數(shù)據(jù)產(chǎn)品。

1.2.2 樣本選擇 本研究結(jié)合研究區(qū)2018年谷歌地球?qū)I(yè)版(Google Earth Pro)高分辨率歷史影像，根據(jù)研究區(qū)內(nèi)主要的地物類型，通過目視解譯的方法，選取研究區(qū)冬小麥—夏玉米、辣椒、大棚蔬菜、其他地物(建設用地、林地等未識別地類)4類樣本，并保證在選擇過程中樣本的隨機性和均勻分布。選擇的樣本數(shù)據(jù)主要分為2個部分用途，一部分用于作物分類方法的時間序列構(gòu)建，一部分用于作物分類結(jié)果的精度驗證，具體樣本選擇見表2。

2 研究方法

2.1 決策樹模型

遙感領域的決策樹分類模型是一種自上而下且根據(jù)不同地物的特殊屬性，結(jié)合先驗知識以及數(shù)據(jù)的統(tǒng)計分析基于訓練集構(gòu)造的1個以二叉樹和多叉樹形式展現(xiàn)的分類規(guī)則模型[21]。決策樹分類模型具有分類規(guī)則靈活、分類數(shù)據(jù)多源、運算高效等特點，在作物識別方面具有廣泛應用[22-24]。本研究利用增強植被指數(shù)(EVI)、比值植被指數(shù)(RVI)、紅光波段反射率、近紅外波段反射率等多源數(shù)據(jù)，共同參與決策樹分類規(guī)則的制定，以更準確地提取不同地物。

表2 樣本個數(shù)統(tǒng)計

2.1.1 MODIS時間序列的構(gòu)建 劉俊偉等認為，在提取多種作物時，尤其是作物周期存在一定重疊的，僅采用單一特征變量來提取作物存在一定困難，可能會出現(xiàn)較多錯分和漏分現(xiàn)象，提取精度會降低，因而通常會選擇多個特征變量作為約束條件進行作物提取[25-26]。而作物的光譜特征和光譜指數(shù)作為表征作物生長過程中重要的指標，可以反映不同作物對波段信號反射的差異性以及作物的生長狀況，故常被組合共同用于作物的提取。因此，本研究選擇EVI、RVI、紅光波段反射率、近紅外波段反射率4個特征變量作為約束條件來提取作物，并基于MODIS時序數(shù)據(jù)利用樣本分別構(gòu)建EVI、RVI、ρR、ρNIR的時序曲線，根據(jù)統(tǒng)計的均值數(shù)據(jù)繪制其變化曲線(圖2)，其中RVI由ρR、ρNIR計算得出。各個時相所對應的日期見表3。

2.1.2 基于MODIS數(shù)據(jù)決策樹模型構(gòu)建 決策樹技術(shù)路線見圖3。結(jié)合表1、圖2可知，區(qū)域內(nèi)對作物識別能夠形成干擾的其他地物主要是建筑、道路等非農(nóng)業(yè)用地，故先將非農(nóng)業(yè)用地從識別影像中剔除掉，觀察發(fā)現(xiàn)非農(nóng)業(yè)用地的EVI值在整個時間序列中都處于較低的水平，且在第7時相非農(nóng)業(yè)用地的EVI值都小于其他作物的EVI值，經(jīng)試驗設置EVI(7)<0.21作為約束條件之一，可以提取出大部分非農(nóng)用地。為保證非農(nóng)業(yè)用地充分剔除掉，引入近紅外波段反射率作為另一個約束條件，由圖2-d 可知，非農(nóng)業(yè)用地在近紅外波段反射率第1時相其反射率遠低于其他地物，且在第1時相非農(nóng)業(yè)用地的最大值都小于其他作物的最小值，因而根據(jù)樣本統(tǒng)計的非農(nóng)業(yè)用地第1時相反射率最大值0.19作為閾值，形成約束條件NIR(1)<0.19。因此，利用約束條件EVI(7)<0.21、NIR(1)<0.19剔除影像內(nèi)的非農(nóng)業(yè)用地。

表3 MODIS數(shù)據(jù)時相

由圖2-a可知，大棚蔬菜在7月下旬，即第14時相其EVI值與其他作物有明顯區(qū)別，這是由于其他作物在7月下旬左右正處于生長旺盛期EVI值較高，因而可以利用該時相設置閾值提取大棚蔬菜，根據(jù)樣本統(tǒng)計值設置約束條件為EVI(14)<0.4。另外，由圖2-b可知，大棚蔬菜的RVI在第14時相其值也遠小于其他作物，為避免提取出來的大棚蔬菜可能混有部分其他作物，將RVI作為另一個約束條件進一步提取大棚蔬菜，經(jīng)反復驗證設置約束條件RVI(14)≤5較合適。因此，最終利用規(guī)則EVI(14)<0.4、RVI(14)≤5來提取大棚蔬菜。

辣椒與冬小麥、夏玉米等作物存在周期重疊問題，單依靠EVI很難完全把辣椒提取出來，而辣椒在紅光波段比其他作物有較強的反射率，因而本研究利用EVI和紅光波段反射率提取該作物。由圖2-a 可知，辣椒在第7時相與冬小麥—夏玉米的EVI值有明顯差距，且辣椒大部分EVI值都集中在0.38附近，可以通過設置合理的EVI區(qū)間閾值將絕大部分辣椒提取出來，根據(jù)樣本統(tǒng)計值設置區(qū)間0.35≤EVI(7)≤0.45提取辣椒較合適。另外，由圖2-c可知，辣椒的紅光波段反射率在春季時較冬小麥、夏玉米都高，故可通過該時期將辣椒與冬小麥—夏玉米區(qū)分開，經(jīng)過試驗發(fā)現(xiàn)設置RED(2)≥0.12作為約束條件可以進一步提取辣椒。因此，選擇提取辣椒的約束條件為0.35≤EVI(7)≤0.45、RED(2)≥0.12。

冬小麥和夏玉米的種植模式為輪種，故提取其中1種作物就可以獲得另一種作物的種植結(jié)構(gòu)。4月上旬是冬小麥生長最旺盛的時期，該時段冬小麥的EVI值明顯高于其他地物，使用第7時相冬小麥EVI最小值0.42作為閾值，大于該閾值的即為冬小麥和少許辣椒。因提取物里面還含有部分辣椒，需要進一步進行提取，因辣椒與冬小麥在第7時相紅光波段反射率有明顯差距，其最小值大于冬小麥的最大值，故使用辣椒紅光波段反射率最小值0.08將兩者分離，即約束條件RED(7)≤0.08。因此，選擇EVI(7)≥0.42、RED(7)≤0.08作為提取冬小麥—夏玉米的條件。

2.2 廣義DEM模型

2.2.1 廣義DEM基本理論 數(shù)字高程模型(DEM)是地形高程數(shù)據(jù)對地形曲面的數(shù)字化表達。通過DEM可以獲取區(qū)域地形因子(坡度、坡向、坡向變率、地形起伏度等)和特征地形要素(山頂點、脊點、山脊線、山谷線等)，使其常應用于土壤、水文、地質(zhì)災害、地貌等多個領域。DEM還可以擴展為廣義DEM，即在二維地理空間上具有連續(xù)變化特征的各種非高程屬性，按照DEM的構(gòu)建方法，使其具備類似地形上的高程屬性，進而可以通過地形分析反映其空間分布特征。

根據(jù)廣義DEM思想[27-28]，可以將遙感植被指數(shù)作為代表植被信息的DEM，進而以植被指數(shù)為基礎計算植被指數(shù)的“地形因子”和“特征地形要素”，根據(jù)這些“地形因子”和“特征地形要素”攜帶的空間信息，可以對農(nóng)作物進行識別和提取。

2.2.2 廣義DEM模型原理和實施流程 選擇提取作物植被指數(shù)最大和最小時期的遙感數(shù)據(jù)做差值運算，通過差值運算增強植被信息，使得不同的地物表現(xiàn)出不同的差值大小，而同一地物的差值將會集中在某一數(shù)值附近，根據(jù)作物的物候信息，作物植被指數(shù)差值將會集中在一個高值區(qū)域，故可選定合適的閾值實現(xiàn)對作物的提取。

廣義DEM模型提取作物基本操作流程見圖4。以冬小麥作物提取為例，其具體實施流程如下：第一，獲取研究區(qū)遙感影像，差值增強植被信息。選擇冬小麥植被指數(shù)最大和最小的兩期影像，并用最大和最小植被指數(shù)相鄰的2個時相做最大值合成和最小值合成，合成后通過差值運算增強農(nóng)作物信息，獲得作物差值增強影像。第二，根據(jù)廣義DEM思想，提取作物差值影像上特征要素“山頂點”作為提取作物的候選樣本。其中“山頂點”提取步驟包括：利用鄰域統(tǒng)計最大值函數(shù)統(tǒng)計差值影像指定鄰域范圍內(nèi)所有輸入像元值的最大值；利用鄰域最大值影像與原始差值影像做差得到“山頂點”影像(山頂點即差值為0的柵格)；利用“山頂點”影像和原始差值影像提取出“山頂點”對應的植被指數(shù)值。第三，“山頂點”提取的作物樣本可能存在一定數(shù)量的錯誤樣本，須要對其進行二次篩選。因而，利用廣義DEM思想中“山頂點”一定位于“山脊線”上的特點，將其作為二次篩選樣本條件，將錯誤樣本去掉。其中，提取山脊線的原理是根據(jù)地面的坡向變率及地面正負地形，取正地形上的坡向變率值較大的部分即為山脊。而坡向變率值較大部分的提取則根據(jù)獲取的坡向變率統(tǒng)計值，選用標準方差作為分類間斷點將坡向變率分為若干類，取最大的一類即為坡向變率值較大的部分。故最終將“山脊線”與“山頂點”疊加進行樣本篩選，取“山頂點”對應的植被指數(shù)的最小值作為作物的提取閾值。

2.2.3 基于MODIS數(shù)據(jù)的廣義DEM方法作物提取

2.2.3.1 數(shù)據(jù)選取與增強 由圖2-a可知，夏玉米在7月下旬生長最旺盛，EVI達到峰值，在其播種期6月和成熟期9月EVI值最小，但由于9月綠色植被生長旺盛容易對夏玉米的提取產(chǎn)生干擾，故本研究選擇7月下旬影像作為最大值、6月中旬影像作為最小值，為增強圖像，使用MODIS第14、15時相進行最大值合成，第11、12時相進行最小值合成，二者差值運算獲得差值增強影像。

辣椒在5月上旬到達生長旺盛期，EVI值達到最大，故可以選擇MODIS中相對應的第8時相作為最大值影像；春植辣椒一般在2月進行播種，此時對應的EVI值最小，所以可以選擇MODIS中相對應的第4時相作為最小值影像。根據(jù)圖像增強處理，選擇MODIS第7、8時相進行最大值合成，第3、4時相進行最小值合成，并進行差值運算獲得增強影像。

由圖2-a還可知，在10月中旬各種作物都基本已經(jīng)成熟或剛播種，所以EVI值較低，而大棚蔬菜在10月中旬EVI值較其他作物都較高，故選擇10月中旬作為大棚蔬菜的最大值影像；觀察EVI時序曲線在10月以后，冬小麥、辣椒的EVI值開始逐漸升高，而大棚蔬菜的EVI逐漸開始降低，在12月左右其他作物EVI值比大棚蔬菜都高，故選擇12月中旬作為大棚蔬菜的最小值影像。因此，選擇MODIS第18、19時相進行最大值合成，第22、23時相進行最小值合成，對圖像做增強處理，進而通過差值運算獲得差值影像。

2.2.3.2 EVI廣義DEM“山頂點”提取 以冬小麥—夏玉米提取為例，根據(jù)廣義DEM“山頂點”提取原理，首先應用鄰域最大值統(tǒng)計算法，采用3×3窗口對冬小麥—夏玉米的差值影像取鄰域最大值獲得最大值影像，用該影像減去原始差值影像，差值結(jié)果像元值為0的即為冬小麥—夏玉米EVI差值影像的“山頂點”，進而提取對應“山頂點”的EVI差值數(shù)值，構(gòu)建冬小麥—夏玉米提取樣本數(shù)據(jù)集。

2.2.3.3 二次篩選樣本點 利用廣義DEM思想，“山頂點”必定位于“山脊線”上，采用3×3的移動窗口計算“坡向”和“坡向變率”，統(tǒng)計“坡向變率”的標準方差，將其作為分類間斷點，提取“坡向變率”值較大部分即為山脊線，將“山脊線”與“山頂點”相交，二次篩選樣本的EVI差值數(shù)值，取其最小值作為冬小麥—夏玉米的提取閾值，其樣本篩選見表4。由于辣椒和冬小麥—夏玉米生長期的重疊，故在提取辣椒時應將冬小麥—夏玉米提取結(jié)果先從提取辣椒影像中掩膜去掉，再進行辣椒和大棚蔬菜的提取，進而對作物辣椒和大棚蔬菜按照相同方法和流程執(zhí)行“2.2.3.2”節(jié)、“2.2.3.3”節(jié)部分的步驟，得到樣本篩選表見表4。

表4 主要作物樣本篩選

3 結(jié)果與分析

3.1 主要作物提取結(jié)果

本研究利用決策樹和廣義DEM分類方法，獲取漳滏河灌區(qū)主要作物種植結(jié)構(gòu)，結(jié)果見圖5、圖6。由圖5和圖6可知，冬小麥—夏玉米是漳滏河灌區(qū)主要的種植作物，其分布范圍廣、種植面積大，且種植區(qū)域基本固定，而辣椒、大棚蔬菜種植面積相對較小且較分散，但這2種作物主要的種植地方都較集中。從決策樹分類結(jié)果來看，辣椒主要位于漳滏河灌區(qū)的東北地區(qū)，即雞澤縣、邱縣和曲周縣，其余地方零散分布；而廣義DEM辣椒分類結(jié)果同樣主要集中在這3縣范圍內(nèi)，但其大部分都集中在雞澤縣附近。賈令晨認為，辣椒是雞澤縣的特色優(yōu)勢產(chǎn)業(yè)和支柱產(chǎn)業(yè)，全縣常年大規(guī)模種植辣椒，已經(jīng)形成了集辣椒種植、加工、銷售完整的產(chǎn)業(yè)體系，因而雞澤縣也被譽為“中國辣椒之鄉(xiāng)”[29]。近幾年，隨著辣椒產(chǎn)業(yè)的發(fā)展也帶動了周邊曲周縣、邱縣、永年區(qū)等地區(qū)開始種植辣椒，逐漸形成了以雞澤縣為中心，其他地區(qū)為輔的辣椒種植基地，這與本研究決策樹和廣義DEM的分類結(jié)果基本相似。

由圖6可知，2個分類結(jié)果大棚蔬菜都主要分布在永年區(qū)和曲周縣，但二者在面積和分布方面存在一定差異，其中決策樹提取的大棚蔬菜都主要集中在永年區(qū)中部地區(qū)，在曲周縣只有少部分，而廣義DEM提取的大棚蔬菜在永年區(qū)中部地區(qū)較少，部分分布于永年區(qū)與曲周縣東北交界處，絕大部分分布在曲周縣中部偏西地區(qū)。

3.2 精度對比

對決策樹和廣義DEM方法的分類結(jié)果精度驗證進行統(tǒng)計，結(jié)果見表5、表6。由表5、表6可知，決策樹的總體分類精度為83.37%，κ系數(shù)為0.78，廣義DEM的總體分類精度為85.33%，κ系數(shù)為0.80，2種方法的總體精度都達到80%以上，分類效果都較理想。但與決策樹分類方法相比，廣義DEM方法的總體分類精度略高于決策樹，較決策樹精度提升了1.96百分點，且各作物的分類精度也都有小幅度提升。因而就提取精度和使用數(shù)據(jù)量的多少而言，廣義DEM在數(shù)據(jù)量較小的情況下也能夠?qū)崿F(xiàn)對作物的識別和區(qū)分。

3.3 方法對比評價

決策樹分類方法在作物識別領域已經(jīng)得到廣泛的驗證和應用，該方法對于各種尺度的范圍、各種時間和空間尺度的數(shù)據(jù)都具有良好的適用性。比較決策樹和廣義DEM分類方法二者各有其特點和優(yōu)勢，決策樹方法對數(shù)據(jù)量需求較大，需要人為采集大量的訓練樣本，且實施過程需要人為干預選取閾值，決策樹分類方法可以將多種數(shù)據(jù)結(jié)合到一起，共同參與決策。相較于決策樹分類方法，廣義DEM分類方法對數(shù)據(jù)量需求較小，人工干預較少，可以自動選取樣本，自動化程度較高，但其通常只能選取1種數(shù)據(jù)作為數(shù)據(jù)源，對于某些生長周期存在重疊、植被指數(shù)數(shù)值相近的作物，提取可能較麻煩。

4 結(jié)論與討論

本研究利用MODIS時間序列數(shù)據(jù)采用決策樹和廣義DEM方法提取漳滏河灌區(qū)的主要作物，并對2種方法進行對比評價，得到以下結(jié)論。第一，廣義DEM作為一種新的作物提取方法，具有很好的可性行和適用性，相對于決策樹而言，不僅對數(shù)據(jù)量要求較小，且人工干預較少，自動化程度較高，有較好的提取效果。第二，廣義DEM方法與決策樹對漳滏河灌區(qū)的作物提取結(jié)果較一致，其中大宗作物冬小麥—夏玉米基本一致，辣椒和大棚蔬菜在面積和分布上存在一定差距，但基本與作物的主要種植區(qū)域吻合。

本研究僅采用高時間、低空間分辨率的MODIS數(shù)據(jù)來探討廣義DEM的可行性和提取優(yōu)勢，未采用中高空間分辨率遙感數(shù)據(jù)對其進行進一步驗證。由于廣義DEM對數(shù)據(jù)量需求較小，在下一步的研究中可利用中高空間分辨率數(shù)據(jù)進行作物提取，以進一步驗證廣義DEM的可行性、適用性以及其提取優(yōu)勢。

表5 決策樹主要作物提取結(jié)果混淆矩陣

表6 廣義DEM方法作物提取結(jié)果混淆矩陣