劉 瑩，朱秀芳，徐 昆

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

0 引 言

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

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

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

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

1 研究區(qū)概況

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

圖1 2017年美國(guó)內(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ù)說(shuō)明

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

1）氣象數(shù)據(jù)為Prism降水?dāng)?shù)據(jù)，來(lái)源于俄勒岡州立大學(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的增強(qiáng)型植被指數(shù)（Enhanced Vegetation Index, EVI）和標(biāo)準(zhǔn)化差值植被指數(shù)（NDVI），MOD16A2實(shí)際蒸散發(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，時(shí)間分辨率為8 d或者16 d。

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

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

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

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

2.2 隨機(jī)森林特征變量評(píng)價(jià)原理

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

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

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

2.3 研究方法

2.3.1 樣本點(diǎn)選擇

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

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

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

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

2.3.2 特征變量計(jì)算

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

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

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

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

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

除了單一元素的特征變量外，本研究還加入了四個(gè)結(jié)合多個(gè)特征運(yùn)算得到的綜合特征變量，包括灌溉概率指數(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 特征變量評(píng)價(jià)和優(yōu)選

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

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

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

3 結(jié)果與分析

3.1 隨機(jī)森林特征評(píng)價(jià)最優(yōu)參數(shù)的取值

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

圖3 分裂屬性個(gè)數(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í)繪制得到的決策樹數(shù)量與袋外樣本誤判率均值的關(guān)系圖如圖4所示。當(dāng)決策樹數(shù)量低于300時(shí)，決策樹數(shù)量的變化會(huì)對(duì)模型的分類效果產(chǎn)生比較顯著的影響；決策樹數(shù)量達(dá)到300后，袋外樣本誤判率均值逐漸趨于穩(wěn)定；當(dāng)決策樹數(shù)量達(dá)到500后，決策樹數(shù)量每增加100，袋外樣本誤判率均值的變化低于0.1個(gè)百分點(diǎn)。因此，決策樹數(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對(duì)備選的82個(gè)特征變量進(jìn)行了重要性評(píng)估，按特征重要性排序依次加入特征進(jìn)行分類，利用袋外樣本進(jìn)行精度評(píng)價(jià)得到的結(jié)果如圖5所示。由圖可知，重要性前4的特征變量（CWSI，EVI，CWSI和IPI）對(duì)提升整體分類精度的影響比較大，重要性5～16的特征變量會(huì)提高模型的分類精度，但是單個(gè)變量的加入對(duì)提升模型分類精度的影響比前幾個(gè)變量有所降低，在加入重要性前16的特征變量后分類結(jié)果基本達(dá)到穩(wěn)定，總體分類精度達(dá)到88.44%。分類的特征變量數(shù)增加到16以后，隨機(jī)森林模型的分類精度呈現(xiàn)小幅度上下波動(dòng)，表示重要性程度排在16位以后的特征變量對(duì)提升模型的分類效果影響不大。因此本研究選取重要性程度前16位的特征變量為最優(yōu)特征變量集合。具體每個(gè)變量的重要性程度及加入變量后的總體分類精度如表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個(gè)最佳特征變量中，包含8個(gè)綜合特征變量、7個(gè)植被特征變量以及1個(gè)土壤特征變量，氣象特征變量沒(méi)有入選。農(nóng)業(yè)干旱指數(shù)（作物水分脅迫指數(shù)）、灌溉概率指數(shù)、植被指數(shù)以及晝夜溫差是區(qū)分灌溉農(nóng)田和雨養(yǎng)農(nóng)田的敏感特征變量。除了灌溉概率指數(shù)不是具備多個(gè)時(shí)相的特征變量外，其他特征變量的不同時(shí)相對(duì)于灌溉農(nóng)田和雨養(yǎng)農(nóng)田的區(qū)分也具有不同的敏感性。對(duì)于作物水分脅迫指數(shù)來(lái)說(shuō)，幾乎生長(zhǎng)季內(nèi)各個(gè)月份以及整個(gè)生長(zhǎng)季的作物水分脅迫指數(shù)（CWSI、CWSI、CWSI、CWSI、CWSI、CWSI和CWSI）都具有較高的敏感性；對(duì)于植被指數(shù)（EVI、EVI、EVI、EVI、NDVI、NDVI和NDVI）來(lái)說(shuō)，對(duì)灌溉農(nóng)田和雨養(yǎng)農(nóng)田的區(qū)分較為敏感的時(shí)相集中在生長(zhǎng)季的后期；對(duì)于晝夜溫差來(lái)說(shuō)，9月份是區(qū)分灌溉農(nóng)田和雨養(yǎng)農(nóng)田的最敏感時(shí)相。

4 討 論

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

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

研究結(jié)果還顯示同一特征變量在不同時(shí)間對(duì)于灌溉的識(shí)別能力不同，整體來(lái)說(shuō)，各個(gè)月份的作物水分脅迫指數(shù)CWSI對(duì)于灌溉的識(shí)別都很有幫助，這也是容易理解的，因?yàn)楣喔茸钪苯拥淖饔镁褪蔷徑廪r(nóng)業(yè)干旱。CWSI由1-AET/PET計(jì)算得到，目前MODIS數(shù)據(jù)提供了AET和PET產(chǎn)品，該數(shù)據(jù)在某些地區(qū)存在缺失，提高該數(shù)據(jù)產(chǎn)品的質(zhì)量有益于推動(dòng)基于該數(shù)據(jù)產(chǎn)品的灌溉耕地產(chǎn)品的生產(chǎn)。

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

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

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

5 結(jié) 論

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

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

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

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