用于灌溉耕地制圖的特征變量優(yōu)選

劉 瑩，朱秀芳，徐 昆

（1. 北京師范大學(xué)遙感科學(xué)國家重點實驗室，北京 100875；2. 北京師范大學(xué)環(huán)境演變與自然災(zāi)害教育部重點實驗室，北京 100875；3. 北京師范大學(xué)地理科學(xué)學(xué)部遙感科學(xué)與工程研究院，北京 100875；4. 山東黃河河務(wù)局山東黃河信息中心，濟南 250013）

0 引 言

灌溉是一種重要的土地管理方式，是用來增加作物產(chǎn)量或者減輕干旱等極端氣候所造成的消極影響的重要手段。在過去的50年里，全球灌溉面積增加了接近一倍。灌溉耕地占全球耕地面積的20%，貢獻了40%的糧食產(chǎn)量。灌溉一方面能夠維持作物產(chǎn)量的穩(wěn)定，促進作物產(chǎn)量的增長，保證面對氣候變化以及干旱、洪水和熱浪等極端氣候事件時的糧食安全，另一方面又會增加土壤濕度和地表蒸散，降低近地面地表溫度，改變熱通量的分割和水分的再分配，對水循環(huán)、能量循環(huán)、生物圈和大氣圈的相互作用以及氣候變化產(chǎn)生間接影響。

明確灌溉的空間范圍、類型以及時序變化有助于開展糧食安全、水資源管理和氣候變化等相關(guān)研究，灌溉耕地制圖的發(fā)展也會為相關(guān)研究提供更多的數(shù)據(jù)基礎(chǔ)和實踐的可能性。目前灌溉耕地制圖的方法可以歸納為兩個類型：基于遙感分類的灌溉耕地制圖和基于統(tǒng)計數(shù)據(jù)空間化的灌溉耕地制圖?；谶b感分類的灌溉耕地制圖通常會利用單個或者一系列光譜特征和變量作為灌溉耕地和其他地物類別的區(qū)分特征，采用監(jiān)督或者非監(jiān)督的分類方法來提取灌溉耕地?；诮y(tǒng)計數(shù)據(jù)空間化的灌溉耕地制圖通常分為兩個步驟進行：1）建立空間分配規(guī)則；2）依據(jù)空間分配規(guī)則，將表征灌溉耕地信息的統(tǒng)計數(shù)據(jù)分配到格網(wǎng)上。第一步中用來建立分配規(guī)則的特征通常是能夠表示灌溉耕地信息相對數(shù)量的特征量，例如灌溉設(shè)施分布、作物復(fù)種指數(shù)等。另外，在統(tǒng)計數(shù)據(jù)空間化方法中，還有一些學(xué)者自定義了一些參數(shù)來進行灌溉范圍的確定。例如，2014年Zhu等利用歸一化植被指數(shù)（Normalized Difference Vegetation Index, NDVI）和降水數(shù)據(jù)依據(jù)長期作物產(chǎn)量高且穩(wěn)定的區(qū)域有更大幾率被灌溉以及灌溉比雨養(yǎng)的作物產(chǎn)量更加穩(wěn)定的假設(shè)創(chuàng)建了三個灌溉潛力指數(shù)，依據(jù)灌溉潛力指數(shù)的大小對灌溉普查數(shù)據(jù)進行空間分配得到2000年中國1 km空間分辨率的灌溉耕地和雨養(yǎng)耕地密度產(chǎn)品。2019年Zohaib等使用衛(wèi)星數(shù)據(jù)集與地表模型模擬的再分析數(shù)據(jù)集之間土壤濕度的正偏差以及地表溫度和地表反照率的負偏差識別灌溉的空間范圍，第一次采用自下而上的方法監(jiān)測了全球?qū)嶋H灌溉區(qū)域。

無論是基于遙感分類的灌溉耕地制圖還是基于統(tǒng)計數(shù)據(jù)空間化的灌溉耕地制圖，構(gòu)建和選擇表征灌溉耕地信息的特征變量都是最重要的環(huán)節(jié)之一。以往大量研究均指出特征變量的選擇對于分類的精度影響很大。例如，朱秀芳等分析了不同分類特征對冬小麥面積測量精度影響，指出紋理特征和植被指數(shù)的加入，不一定能提高分類精度，最優(yōu)特征和研究區(qū)本身的特點以及分類器等都有關(guān)系。研究者往往在分類前會對分類特征進行優(yōu)選。例如，劉杰等基于多時相Landsat 8 OLI數(shù)據(jù)，提取光譜、紋理、植被指數(shù)等217個特征波段，利用隨機森林中的基尼系數(shù)（GINI Coefficient）優(yōu)選特征子集，在新疆維吾爾自治區(qū)阿克蘇地區(qū)溫宿縣進行了作物類型的精細識別。朱秀芳等利用灰度共生矩陣提取無人機影像的紋理特征，使用隨機森林算法進行紋理特征優(yōu)選，在此基礎(chǔ)上提取了覆膜農(nóng)田的面積和分布。王庚澤等提出了改進分離閾值組合式特征優(yōu)選算法，從多時相的Sentinel 2原始光譜波段、不同時相波段差值和比值、傳統(tǒng)植被指數(shù)和紅邊植被指數(shù)、第一主成分的紋理特征合計183個特征中選擇最佳特征，進而利用隨機森林算法提取河北省藁城區(qū)的秋糧作物。

目前常用的灌溉耕地制圖的特征變量中，各變量對灌溉耕地識別精度的貢獻并不相同。明確變量貢獻大小可以指導(dǎo)后續(xù)其他地區(qū)進行灌溉耕地制圖特征變量的選擇，進而提高灌溉耕地的特征變量的選擇效率和最終的識別精度。為此，本研究基于隨機森林，在美國內(nèi)布拉斯加州，對灌溉耕地制圖中4類82個特征變量進行比較分析，優(yōu)選最佳特征變量，為后續(xù)相關(guān)研究中特征變量的選擇提供參考。

1 研究區(qū)概況

本研究選擇美國的內(nèi)布拉斯加州（Nebraska）作為研究區(qū)（圖1），該州緯度位置位于40°N～43°N之間，經(jīng)度位置位于95°W～105°W之間。內(nèi)布拉斯加州有39%的土地面積被用于農(nóng)業(yè)生產(chǎn)，是美國主要的農(nóng)業(yè)生產(chǎn)地區(qū)之一。該州有超過10萬的活躍灌溉水井，并且每十年的灌溉水井增量達1萬個，是世界上灌溉最密集的區(qū)域之一。該州主要的種植作物包括玉米（2017年收獲面積占比50%）和大豆（2017年收獲面積占比30%），其他次要作物包括冬小麥，高粱和苜蓿（2017年總收獲面積占比約11%）。4－10月是該州主要作物的生長季。近年來有學(xué)者制作了包含該地區(qū)的長時間序列中高分辨率灌溉耕地分布圖。選擇該州作為研究區(qū)便于利用已有的灌溉耕地數(shù)據(jù)產(chǎn)品選擇樣本數(shù)據(jù)進行灌溉特征變量的對比分析。

圖1 2017年美國內(nèi)布拉斯加州灌溉農(nóng)田面積占農(nóng)田面積比例Fig.1 Proportion of irrigated farmland area to farmland area in Nebraska, United States in 2017

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)說明

本研究的輸入數(shù)據(jù)主要包括4個部分：氣象數(shù)據(jù)、遙感數(shù)據(jù)、灌溉數(shù)據(jù)和輔助數(shù)據(jù)（表1）。數(shù)據(jù)的基本信息和預(yù)處理過程如下：

1）氣象數(shù)據(jù)為Prism降水數(shù)據(jù)，來源于俄勒岡州立大學(xué)的 PRISM 氣候研究組，下載網(wǎng)址為：http://prism.oregonstate.edu/recent/。該數(shù)據(jù)的空間分辨率為0.04°，空間參考為GCS North American 1983。

2）遙感數(shù)據(jù)包括MOD09A1 反射率，MOD11A2 地表溫度數(shù)據(jù)（Land Surface Temperature, LST），MOD13A1的增強型植被指數(shù)（Enhanced Vegetation Index, EVI）和標準化差值植被指數(shù)（NDVI），MOD16A2實際蒸散發(fā)（Actual Evapotranspiration, AET）和潛在蒸散發(fā)（Potential Evapotranspiration, PET），下載網(wǎng)址為https://lpdaac.usgs.gov/products/。這些數(shù)據(jù)的空間參考類型為GCS WGS 1984，空間分辨率為500 m或者1 000 m，時間分辨率為8 d或者16 d。

3）灌溉數(shù)據(jù)包括美國高平原含水層年度灌溉圖（Annual Irrigation Maps-High Plains Aquifer，AIM-HPA）、中分辨率成像光譜儀灌溉農(nóng)業(yè)數(shù)據(jù)集（The Moderate Resolution Imaging Spectroradiometer Irrigated Agriculture Dataset for the United States，MIrAD-US）、內(nèi)布拉斯加州灌溉設(shè)施分布圖和已注冊灌溉水井矢量數(shù)據(jù)。AIM-HPA是覆蓋美國高平原含水層的1984－2017年每年30 m空間分辨率的灌溉耕地數(shù)據(jù)集。該數(shù)據(jù)集是由Deines等利用Landsat影像、環(huán)境變量和地面實況數(shù)據(jù)，基于隨機森林分類器分類，并在分類后應(yīng)用貝葉斯土地覆蓋更新算法填補圖像空白，降低灌溉時間序列中的誤差得到的。MIrAD-US是由美國地質(zhì)調(diào)查局地球資源觀測與科學(xué)中心開發(fā)的與美國農(nóng)業(yè)普查數(shù)據(jù)同步更新的美國灌溉農(nóng)業(yè)數(shù)據(jù)集。該數(shù)據(jù)集是利用統(tǒng)計數(shù)據(jù)空間化的方法建模得到，以統(tǒng)計數(shù)據(jù)作為總量約束，以年NDVI峰值作為分配依據(jù)，將灌溉面積統(tǒng)計數(shù)據(jù)分配到空間格網(wǎng)上得到的。該數(shù)據(jù)集包含2002、2007、2012和2017年空間分辨率為250和1 000 m的灌溉耕地數(shù)據(jù)。內(nèi)布拉斯加州灌溉設(shè)施分布圖由內(nèi)布拉斯加大學(xué)林肯分校先進土地管理信息技術(shù)中心提供。該數(shù)據(jù)圖層以面狀矢量文件的形式提供了根據(jù)2005年多時相的Landsat 5衛(wèi)星影像和農(nóng)業(yè)服務(wù)局的正射影像確定的內(nèi)布拉斯加州2005年生長季的中心樞紐灌溉系統(tǒng)和其他灌溉系統(tǒng)分布圖。內(nèi)布拉斯加州已注冊灌溉水井來源于內(nèi)布拉斯加州自然資源部網(wǎng)站提供的注冊地下水井數(shù)據(jù)。該數(shù)據(jù)以矢量點的形式存儲，利用其屬性表篩選出用途為灌溉且狀態(tài)為活躍的水井構(gòu)成已注冊活躍灌溉水井數(shù)據(jù)。

4）輔助數(shù)據(jù)包括農(nóng)作物圖層、作物物候和研究區(qū)矢量。農(nóng)作物圖層是由美國農(nóng)業(yè)部的美國國家農(nóng)業(yè)統(tǒng)計局提供的農(nóng)作物數(shù)據(jù)圖層，空間分辨率為30 m。作物物候數(shù)據(jù)來源于美國農(nóng)業(yè)部的美國國家農(nóng)業(yè)統(tǒng)計局于2010年10月發(fā)布的農(nóng)作物播種和收獲日期。該數(shù)據(jù)主要被用來根據(jù)內(nèi)布拉斯加州主要農(nóng)作物的播種和收獲時間確定研究中待分析數(shù)據(jù)的候選時間：4—10月。研究區(qū)內(nèi)布拉斯加州的矢量邊界數(shù)據(jù)來源于GADM (Database of Global Administrative Areas)。GADM提供所有國家以及地區(qū)的行政邊界圖，其下載網(wǎng)址為http://www.gadm.org/。

表1 數(shù)據(jù)說明Table 1 Data description

將所有氣象數(shù)據(jù)和遙感數(shù)據(jù)的空間分辨率統(tǒng)一到500 m，時間分辨率統(tǒng)一到月尺度，空間參考統(tǒng)一為Albers Conic Equal Area。將其他數(shù)據(jù)的空間參考統(tǒng)一為Albers Conic Equal Area。

2.2 隨機森林特征變量評價原理

隨機森林由Breiman于2001年提出，它是一個包含多個決策樹的集成分類器，其中的每棵決策樹之間都是沒有關(guān)聯(lián)的。對于每個輸入樣本，每棵決策樹都會給出一個分類結(jié)果，最后取所有決策樹中分類結(jié)果最多的類別作為這個輸入樣本的最終分類結(jié)果。隨機森林算法分類精度高、處理速度快、可處理高維數(shù)據(jù)和多重共線性數(shù)據(jù)且對過度擬合不敏感，是遙感分類的常用算法之一，也曾被成功應(yīng)用于灌溉耕地制圖的研究當中。

隨機森林算法還可以對特征變量的重要程度和貢獻進行分析評價，以保證使用最少的特征數(shù)量達到最佳的分類精度，通過特征變量的篩選，可以減少分類數(shù)據(jù)的冗余、降低不可靠變量對分類效果的影響、提高遙感分類的工作效率，并為后續(xù)相關(guān)研究進行特征變量的選擇提供依據(jù)。在隨機森林算法中主要有兩個指標可以對特征進行重要性評估：一種是基于袋外誤差的平均準確度降低值（Mean Decrease Accuracy, MDA）；另一種是基于基尼不純度的平均不純度降低值（Mean Decrease Gini, MDG）。本文選擇與分類精度關(guān)系更加密切的MDA進行特征重要性的評估。隨機森林每棵決策樹的生成過程中都有接近三分之一的樣本沒有被使用，這些樣本被稱為袋外樣本。MDA主要是利用這些袋外樣本來對變量的重要性進行評估，其計算過程大致分為以下3個步驟：1）對于每棵決策樹，袋外樣本數(shù)據(jù)的預(yù)測錯誤率都會被記錄下來；2）將袋外樣本數(shù)據(jù)的特征變量的值隨機打亂（這個步驟相當于依次將待評估的特征變量替換成噪聲），再次記錄下來袋外樣本數(shù)據(jù)的預(yù)測錯誤率；3）對于每棵決策樹，計算這兩次袋外樣本預(yù)測錯誤率的差值，將所有決策樹上的差值求均值。其表達公式如下：

式中ntree代表這個森林中樹的數(shù)量，E代表第棵樹在對X特征變量的值隨機打亂之前的袋外樣本錯誤率，EP代表第棵樹在對X特征變量的值隨機打亂之后的袋外樣本錯誤率。MDA的值越大，表示當前特征被替換成噪聲后對模型預(yù)測精度的改變程度越大，也就是這個特征的重要性越高。

2.3 研究方法

2.3.1 樣本點選擇

結(jié)合作物圖層CDL、已有的灌溉耕地分布圖AIM-HPA和MIrAD、內(nèi)布拉斯加州灌溉設(shè)施分布圖和已注冊灌溉水井分布圖，按照下述規(guī)則在研究區(qū)內(nèi)選取灌溉耕地樣本點和雨養(yǎng)耕地樣本點。

灌溉耕地樣本點的選取規(guī)則分為兩類：第一類適用于AIM-HPA灌溉耕地分布圖有覆蓋的區(qū)域；第二類適用于AIM-HPA沒有覆蓋的內(nèi)布拉斯加州東部區(qū)域（圖2）。第一類選取規(guī)則需要同時滿足以下條件：1）AIM-HPA在2002－2017年均為灌溉耕地；2）MIrAD在2002、2007、2012和2017年均為灌溉耕地；3）2005年具有灌溉設(shè)備的區(qū)域。由于AIM-HPA沒有覆蓋的內(nèi)布拉斯加州東部區(qū)域?qū)儆跐駶櫟貐^(qū)，該地區(qū)的水分供給比較充足，對于灌溉的需求較少，經(jīng)過疊加分析發(fā)現(xiàn)該部分地區(qū)MIrAD在2002、2007、2012和2017年四年均為灌溉耕地的像元也比較少，因此第二類選取更改為需要同時滿足以下條件：1）MIrAD在2017年為灌溉耕地；2）2005年具有灌溉設(shè)備的區(qū)域；3）靠近已注冊活躍灌溉水井點。

圖2 高平原含水層年度灌溉圖（AIM-HPA）及在研究區(qū)中的覆蓋情況Fig.2 Annual Irrigation Maps-High Plains Aquifer (AIM-HPA)and its coverage in the study area

由于AIM-HPA和MIrAD灌溉作物分布圖只包含灌溉作物像元，因此需要結(jié)合作物圖層CDL對雨養(yǎng)耕地樣本進行提取。對于30 m空間分辨率的AIM-HPA，采用重新編碼后的30 m分辨率的作物圖層進行雨養(yǎng)耕地的提?。簩IM-HPA為非灌溉耕地像元且CDL為作物的像元識別為雨養(yǎng)耕地。對于250 m空間分辨率的MIrAD，采用重新編碼且進行眾數(shù)聚合后的250 m分辨率的作物圖層進行雨養(yǎng)耕地的提?。簩IrAD為非灌溉耕地像元且CDL為作物的像元識別為雨養(yǎng)耕地。雨養(yǎng)耕地樣本點的選取規(guī)則也按照AIM-HPA有無覆蓋分為兩類。在AIM-HPA有覆蓋的區(qū)域需要同時滿足以下條件：1）AIM-HPA在2002－2017年均為雨養(yǎng)耕地；2）MIrAD在2002、2007、2012和2017年均為雨養(yǎng)耕地；3）2005年不具有灌溉設(shè)備的區(qū)域。在AIM-HPA沒有覆蓋的區(qū)域需要同時滿足以下條件：1）MIrAD在2017年為雨養(yǎng)耕地；2）2005年不具有灌溉設(shè)備的區(qū)域；3）遠離已注冊活躍灌溉水井點。

依據(jù)以上規(guī)則選擇了440個灌溉耕地樣本點和343個雨養(yǎng)耕地樣本點。

2.3.2 特征變量計算

本研究選取一些常用于表征耕地受灌溉可能性的特征變量作為遙感分類中的候選特征變量，包括氣象特征變量、植被特征變量、土壤特征變量以及一些經(jīng)過多個特征運算得到的綜合特征變量，其中一些變量包含了生長季內(nèi)每個月份以及整個生長季內(nèi)的值，總計共82個特征變量，每個特征變量的計算方法以及其他信息如表2所示。

表2 待分析的特征變量Table 2 Feature variables to be analyzed

氣象特征變量包括能對作物可用水分產(chǎn)生直接影響的降水量（Precipitation, Pre）數(shù)據(jù)。氣象條件是決定農(nóng)作物生長狀態(tài)的基礎(chǔ)條件，也是決定是否需要進行灌溉的先決條件。

植被特征變量主要包含植被的生產(chǎn)力水平以及植被的含水量情況。其中，生產(chǎn)力水平是利用植被指數(shù)來表示的。植被指數(shù)是監(jiān)測作物長勢以及生產(chǎn)力水平的有效手段，也是作物所處的氣候條件以及是否受到灌溉等人為管理措施的間接體現(xiàn)。當氣候條件適宜時，灌溉是增加作物可用水分、提高植被生產(chǎn)力的有效方式；當作物受到干旱等極端氣候條件的脅迫時，灌溉可以通過改變作物產(chǎn)量波動的閾值或者使得作物產(chǎn)量與氣候條件解耦來保持作物生產(chǎn)力水平的穩(wěn)定性。因此，可以通過植被指數(shù)反映的作物生長狀況來間接反映作物的生長過程中是否有灌溉的參與。本研究共選取了三種植被指數(shù)作為候選的植被特征變量：歸一化植被指數(shù)（NDVI）、增強型植被指數(shù)和植被綠度指數(shù)（Greenness Index, GI）。植被的含水率情況是用歸一化水指數(shù)（Normalized Difference Water Index, NDWI）表示的，該指數(shù)可以反映植被冠層的液態(tài)水含量。

土壤特征變量利用地表溫度（Land Surface Temperature，LST）表示。研究表明，灌溉可以提高土壤濕度，增加農(nóng)田的蒸散量，改變地表屬性，這些生物物理影響也會進一步影響到地表能量的分割，反饋到局部氣候上使得地表溫度降低。因此，地表溫度也可以作為間接反映作物是否被灌溉的特征變量。本研究中選取了日地表溫度（LSTday）、夜地表溫度（LSTnight）以及晝夜溫差（LSTdifference）作為分類過程中候選的土壤特征變量。

除了單一元素的特征變量外，本研究還加入了四個結(jié)合多個特征運算得到的綜合特征變量，包括灌溉概率指數(shù)（Irrigation Probability Index, IPI）、作物水分虧缺指數(shù)（Crop Water Deficit Index, CWDI）、作物水分脅迫指數(shù)（Crop Water Stress Index, CWSI）和水分調(diào)整綠度指數(shù)（Water-adjusted green index, WGI）。

2.3.3 特征變量評價和優(yōu)選

隨機森林算法需要進行兩個重要參數(shù)的設(shè)置：決策樹的數(shù)量（ntree）和樹木生長時為達到最佳分割效果選擇的分裂屬性個數(shù)（mtry）。理論以及經(jīng)驗上的研究表明：分類精度對于決策樹數(shù)量的敏感程度低于分裂屬性個數(shù)，由于隨機森林分類器的計算效率較高而且不會過擬合，決策樹的數(shù)量可以設(shè)置的盡可能大。通常，選擇一個使得分類誤差達到穩(wěn)定的決策樹數(shù)量即可。分裂屬性個數(shù)通常會被設(shè)置為輸入變量數(shù)的平方根。

本研究利用R語言的randomforest程序包實現(xiàn)隨機森林分類算法。首先，將決策樹數(shù)量設(shè)置為一個較大的數(shù)值（設(shè)置ntree=1 000），保證在這個數(shù)量下的分類誤差已經(jīng)達到穩(wěn)定，依次設(shè)置分裂屬性mtry=1，2，3，…，81（變量總個數(shù)為82）進行試驗，得到模型的袋外樣本誤判率均值隨著mtry的變化情況，找出模型袋外樣本誤判率均值達到最低時的mtry取值作為mtry最優(yōu)值。固定mtry值為最優(yōu)值，繪制ntree=1，10，20，50，100，200，300，400，500，600，700，800，900，1 000時，決策樹數(shù)量與袋外樣本誤判率均值的關(guān)系圖，選擇袋外樣本誤判率均值的變化低于0.1%的臨界值，將其作為ntree的最優(yōu)取值。

依據(jù)確定的適用于特征選擇過程的最佳mtry和ntree參數(shù)，利用MDA對備選的82個特征變量進行了重要性評估。將此82個變量按照MDA重要性降序排列，從重要性最高的變量開始，依次添加變量進行隨機森林分類（由于該步驟中變量數(shù)在隨時變化，因此將mtry設(shè)置為變量數(shù)的平方根，ntree仍然設(shè)置為特征選擇過程中得到的最佳ntree值），利用袋外樣本的總體分類精度對分類效果進行評價，找出分類精度最高時對應(yīng)的特征變量，作為優(yōu)選出的特征變量集合。

3 結(jié)果與分析

3.1 隨機森林特征評價最優(yōu)參數(shù)的取值

當決策樹數(shù)量固定為1 000時，模型的袋外樣本誤判率均值隨著mtry的變化情況結(jié)果如圖3所示。mtry的取值對模型誤判率的影響整體波動不超過4%，當mtry=11時，模型袋外樣本誤判率均值達到最低，為11.33%，當mtry=9，與輸入變量數(shù)的平方根接近時，模型袋外樣本誤判率均值次低，為11.36%。因此，分裂屬性個數(shù)mtry的最優(yōu)取值為11。

圖3 分裂屬性個數(shù)（mtry）與袋外樣本誤判率均值的關(guān)系Fig.3 The relationship between the number of split attributes(mtry) and the mean out-of-bag error

mtry參數(shù)固定為11，ntree為1、10、20、50、100、200、300、400、500、600、700、800、900、1 000時繪制得到的決策樹數(shù)量與袋外樣本誤判率均值的關(guān)系圖如圖4所示。當決策樹數(shù)量低于300時，決策樹數(shù)量的變化會對模型的分類效果產(chǎn)生比較顯著的影響；決策樹數(shù)量達到300后，袋外樣本誤判率均值逐漸趨于穩(wěn)定；當決策樹數(shù)量達到500后，決策樹數(shù)量每增加100，袋外樣本誤判率均值的變化低于0.1個百分點。因此，決策樹數(shù)量ntree的最優(yōu)取值為500。

圖4 決策樹數(shù)量（ntree）與袋外樣本誤判率均值的關(guān)系Fig.4 The relationship between the number of decision trees(ntree) and the mean out-of-bag error

3.2 分類特征優(yōu)選結(jié)果

設(shè)定mtry=11，ntree=500，利用MDA對備選的82個特征變量進行了重要性評估，按特征重要性排序依次加入特征進行分類，利用袋外樣本進行精度評價得到的結(jié)果如圖5所示。由圖可知，重要性前4的特征變量（CWSI，EVI，CWSI和IPI）對提升整體分類精度的影響比較大，重要性5～16的特征變量會提高模型的分類精度，但是單個變量的加入對提升模型分類精度的影響比前幾個變量有所降低，在加入重要性前16的特征變量后分類結(jié)果基本達到穩(wěn)定，總體分類精度達到88.44%。分類的特征變量數(shù)增加到16以后，隨機森林模型的分類精度呈現(xiàn)小幅度上下波動，表示重要性程度排在16位以后的特征變量對提升模型的分類效果影響不大。因此本研究選取重要性程度前16位的特征變量為最優(yōu)特征變量集合。具體每個變量的重要性程度及加入變量后的總體分類精度如表3所示。

圖5 特征數(shù)量與分類精度的關(guān)系Fig.5 Relationship between the number of features and classification accuracy

表3 特征重要性排序及分類總體精度Table 3 Feature importance ranking and overall accuracy of classification

在選取出的16個最佳特征變量中，包含8個綜合特征變量、7個植被特征變量以及1個土壤特征變量，氣象特征變量沒有入選。農(nóng)業(yè)干旱指數(shù)（作物水分脅迫指數(shù)）、灌溉概率指數(shù)、植被指數(shù)以及晝夜溫差是區(qū)分灌溉農(nóng)田和雨養(yǎng)農(nóng)田的敏感特征變量。除了灌溉概率指數(shù)不是具備多個時相的特征變量外，其他特征變量的不同時相對于灌溉農(nóng)田和雨養(yǎng)農(nóng)田的區(qū)分也具有不同的敏感性。對于作物水分脅迫指數(shù)來說，幾乎生長季內(nèi)各個月份以及整個生長季的作物水分脅迫指數(shù)（CWSI、CWSI、CWSI、CWSI、CWSI、CWSI和CWSI）都具有較高的敏感性；對于植被指數(shù)（EVI、EVI、EVI、EVI、NDVI、NDVI和NDVI）來說，對灌溉農(nóng)田和雨養(yǎng)農(nóng)田的區(qū)分較為敏感的時相集中在生長季的后期；對于晝夜溫差來說，9月份是區(qū)分灌溉農(nóng)田和雨養(yǎng)農(nóng)田的最敏感時相。

4 討 論

本研究的結(jié)果進一步證實了遙感分類特征對于待分類目標的識別精度具有很大的影響，從表3可以看出，使用重要性排在前16的變量進行灌溉農(nóng)田制圖時比僅使用重要性排在前2位的變量進行灌溉農(nóng)田制圖的精度提高了12.49個百分點，加入更多的重要性排序靠后的特征對分類精度的提升沒有作用，甚至起到負作用。因此，在進行灌溉制圖之前有必要進行特征選擇工作。

由特征選擇的結(jié)果可以看出，綜合特征變量對灌溉農(nóng)田的識別最有幫助，而氣象特征變量對灌溉農(nóng)田制圖的幫助最小。灌溉可以緩解作物水分脅迫，減緩或者抑制農(nóng)業(yè)干旱的發(fā)生，進而促進作物生長，相比雨養(yǎng)作物同類的灌溉作物的植被生產(chǎn)力更高，灌溉也會影響地表參數(shù)，例如使地表溫度降低、蒸散增強、土壤濕度增加。缺少降水只能決定灌溉需求，但不能決定是否真的發(fā)生了灌溉，因此降水特征變量在本研究中表現(xiàn)出最低的重要性。而能反映灌溉效果或灌溉可能性的指標，如作物水分脅迫指數(shù)CWSI，灌溉概率指數(shù)IPI和植被指數(shù)更有助于判斷是否發(fā)生了灌溉。簡單來說從灌溉結(jié)果出發(fā)比從灌溉的誘因出發(fā)去選擇特征變量能更有效地判定是否發(fā)生了灌溉。

研究結(jié)果還顯示同一特征變量在不同時間對于灌溉的識別能力不同，整體來說，各個月份的作物水分脅迫指數(shù)CWSI對于灌溉的識別都很有幫助，這也是容易理解的，因為灌溉最直接的作用就是緩解農(nóng)業(yè)干旱。CWSI由1-AET/PET計算得到，目前MODIS數(shù)據(jù)提供了AET和PET產(chǎn)品，該數(shù)據(jù)在某些地區(qū)存在缺失，提高該數(shù)據(jù)產(chǎn)品的質(zhì)量有益于推動基于該數(shù)據(jù)產(chǎn)品的灌溉耕地產(chǎn)品的生產(chǎn)。

本研究針對灌溉農(nóng)田制圖的特征變量選擇問題進行了有益的探索，相關(guān)結(jié)果可以為其他研究者篩選灌溉農(nóng)田制圖的特征變量提供參考。但本研究也存在一些不足之處：

1）本研究的研究區(qū)為美國的內(nèi)布拉斯加州，其主要的種植作物為玉米和大豆，尚未在其他種植結(jié)構(gòu)不同的地區(qū)進行分析。中國是農(nóng)業(yè)大國，其灌溉用水量和灌溉面積均位居全球前列。本文所提出的方法在中國有很大的應(yīng)用潛力，然而中國的耕地地塊相比美國更加破碎，種植結(jié)構(gòu)更加復(fù)雜，且南北方差異大，在干旱和半干旱區(qū)為了節(jié)約用水，往往以補給灌溉為主，相比充分灌溉，補給灌溉的信號更弱，可能會影響識別的效果。未來還需要選擇南北方種植結(jié)構(gòu)不同的典型灌區(qū)，進一步的驗證本文所提出的方法。

2）特征選擇方法大致分為過濾式（Filter）和封裝式（Wrapper）兩種。前者在數(shù)據(jù)預(yù)處理步驟中對特征排序，設(shè)定閾值選擇最優(yōu)特征子集，排序準則有相關(guān)系數(shù)、互信息等。后者將分類特征與學(xué)習(xí)算法相結(jié)合，根據(jù)準確率評價每個特征子集，從而選擇最優(yōu)特征子集。本文所用的隨機森林就是封裝式特征選擇方法的一種。與過濾式特征選擇方法相比，隨機森林對噪聲數(shù)據(jù)和存在缺失值的數(shù)據(jù)具有較好的魯棒性，其預(yù)測能力不受多重共線性影響。但是一些具體多重共線性的特征（如NDVI和EVI）的重要性會被相互抵消，從而影響對特征變量的解釋性。對于高度相關(guān)的特征變量會存在信息冗余。未來考慮綜合使用過濾式和封裝式的特征選擇方法，以減少特征變量的信息冗余和提高特征變量的可解釋能力，找到普適性的灌溉耕地制圖特征變量，服務(wù)于大范圍灌溉耕地制圖產(chǎn)品的生產(chǎn)。

5 結(jié) 論

本文選擇有良好灌溉信息數(shù)據(jù)基礎(chǔ)的美國內(nèi)布拉斯加州為研究區(qū)，基于已有耕地空間分布圖和灌溉信息數(shù)據(jù)，提取灌溉耕地樣本和雨養(yǎng)耕地樣本，計算了樣本的82個特征變量，利用隨機森林對比分析了82個特征變量對灌溉耕地識別的重要性，得到如下主要結(jié)論：

1）5月作物水分脅迫指數(shù)，7月增強型植被指數(shù)，4月作物水分脅迫指數(shù)和灌溉概率指數(shù)是重要性前4的特征變量，對提升灌溉耕地制圖精度的影響最為明顯。利用重要性前16的特征變量分類得到的總體分類精度最高，為88.44%。

2）四類特征變量中，對灌溉耕地識別的貢獻程度由大到小為綜合特征變量、植被特征變量、土壤特征變量、氣象特征變量。

3）不同特征變量識別灌溉農(nóng)田的最佳時相也存在差異。生長季后期的植被指數(shù)和9月的晝夜溫差相比其他月份的植被指數(shù)和晝夜溫差更有利于灌溉耕地識別。而對于農(nóng)業(yè)干旱指數(shù)來說，幾乎生長季內(nèi)所有的月份都有助于灌溉耕地的識別。