朱德海 劉逸銘 馮權瀧,3 歐 聰 郭 浩 劉建濤

(1.中國農業(yè)大學土地科學與技術學院， 北京 100083； 2.自然資源部農用地質量與監(jiān)控重點實驗室， 北京 100193；3.中國農業(yè)大學資源與環(huán)境學院， 北京 100193； 4.山東建筑大學測繪地理信息學院， 濟南 250101)

0 引言

農業(yè)部于1988年提出建設“菜籃子工程”，集約化蔬菜生產(以農業(yè)大棚蔬菜為主要代表)已經成為蔬菜生產的主要發(fā)展方向[1]。山東省是我國蔬菜種植大省，同時也是全國設施種植蔬菜的中心。因此，監(jiān)測山東省農業(yè)大棚的空間分布和動態(tài)變化，可為相關部門的政策制定提供技術支撐，并有效促進“菜籃子工程”的順利實施。

目前，農業(yè)大棚數據的獲取多采用抽樣調查和統(tǒng)計上報等方法，由于其工作流程復雜、人工干預過多等問題，導致數據獲取的準確性和時效性較低。相比于傳統(tǒng)統(tǒng)計方法，遙感具有宏觀、動態(tài)、快速等特點，可以實現大面積動態(tài)同步觀測，已用于農業(yè)大棚的信息提取。其中，目視解譯雖有較高的精度，但費時費力[2-3]，難以完成長時間序列、大空間尺度的農業(yè)大棚監(jiān)測任務。因此，基于機器學習、云計算等前沿技術，實現長時間、大范圍、高精度的農業(yè)大棚信息自動提取，成為亟待解決的重要課題。

在農業(yè)大棚自動提取方面，國外開展研究較早，其中AGUILAR等[4]基于WorldView-2和Landsat-8 OLI影像數據，采用面向對象方法和決策樹分類器在西班牙南部某地區(qū)進行了大棚提取試驗，分類精度達到93%。LANORTE等[5]結合土地利用數據和Landsat-8影像數據，提出一種用于提取大棚的塑料污染指數，其分類精度達到94.54%。FRANCISCO等[6]研究表明，在QuickBird和IKONOS等高分辨率衛(wèi)星影像中添加紋理特征，可以有效提高塑料大棚的分類精度。黃振國等[7]提出了一種基于知識分類的大棚提取方法，以山東省壽光市為研究區(qū)，分類精度達到了92.01%。王志盼等[8]以廣東省江門市作為研究區(qū)，提出了一種基于增強型水體指數的大棚分類方法，分類精度達75%。YANG等[9]以山東省濰坊市西北部為研究區(qū)，通過研究Landsat ETM+各波段對大棚的敏感度，提出了大棚指數，其分類精度達到91.2%。孫鈺等[10]采用赤峰市王爺府鎮(zhèn)地區(qū)的無人機航空影像，基于全卷積神經網絡構建了大棚提取模型，其分類精度達到97%。此外還有其他學者提出了相關的大棚提取方法[11-13]。然而，上述研究主要集中在小范圍的典型區(qū)域，較少涉及到大空間尺度、長時間序列的農業(yè)大棚信息提取。

大空間尺度、長時間序列遙感數據集的存儲、組織和計算對軟硬件環(huán)境提出了很高的要求。針對上述問題，谷歌發(fā)布了Google Earth Engine云平臺，可實現全球尺度地球科學數據(尤其是遙感數據，如Landsat、MODIS、Sentinel等)的在線可視化計算和分析處理，從而為大空間尺度、長時間序列遙感數據的實時處理提供了強有力的技術支撐。相比于ENVI等桌面端遙感圖像處理軟件，GEE提供了基于Web端的JavaScript或Python算法開發(fā)工具，可快速實現批處理和可視化等功能[14-16]。國內外學者利用GEE平臺開展了一系列研究，其中孟夢等[17]基于GEE制作了白洋淀地區(qū)1987—2017年的NDVI、NDWI數據集，并對白洋淀的濕地生態(tài)系統(tǒng)進行了景觀分類。MAHDIANPARI等[18]利用GEE平臺和Sentinel-1/2數據對加拿大紐芬蘭地區(qū)進行了濕地制圖，其制圖精度為88.37%。

本文以山東省為研究區(qū)，通過在GEE云平臺上獲取1990、1998、2009、2018年多時相Landsat影像，分析山東省近30年農業(yè)大棚的空間分布和時空動態(tài)變化。

1 研究區(qū)域與數據來源

1.1 研究區(qū)域概況

本文選取山東省(34°22′～38°24′N,114°19′～122°43′E)作為研究區(qū)，如圖1所示。山東省年平均氣溫為11～14℃，年均降水量為550～950 mm。山東省是全國糧食作物和經濟作物重點產區(qū)，素有“糧棉油之庫，水果水產之鄉(xiāng)”之稱[19]。

圖1 研究區(qū)Fig.1 Map of research area

1.2 主要數據源

Landsat衛(wèi)星是由美國地質調查局和美國宇航局共同實施，從1972年至今持續(xù)對地觀測，保存了大量珍貴的歷史影像資料。Landsat衛(wèi)星以14 d左右的重返周期對地球表面進行30 m分辨率成像，包括多光譜和熱量數據，廣泛應用于土地調查、環(huán)境監(jiān)測等領域[20]。

山東省于20世紀80年代末、90年代初開始大力發(fā)展農業(yè)大棚，且在GEE平臺獲取的1990年山東省Landsat影像質量較好，故本文以1990年作為研究起始時間，共選取4個時期來表征近30年的時空變化情況，分別為1990、1998、2009年Landsat-5 TM影像和2018年Landsat-8 OLI影像，由于在GEE中，山東省2008年影像云量較多、質量較差，故選取2009年的數據。Landsat數據在GEE中提供如大氣上層反射率、地表反射率等即用型產品。本文基于GEE平臺完成長時間序列Landsat-5/8遙感影像的輻射定標、幾何校正、云雪去除、影像拼接等處理。同時，本文利用GEE平臺多時相Landsat影像制作了長時間序列NDVI、MNDWI、SAVI等光譜特征和紋理特征數據集。其他輔助分析數據主要包括地形和坡度數據等。

2 研究方法

2.1 工作流程

工作流程如圖2所示，主要步驟如下：

(1)影像調取及預處理：在GEE云平臺上分別獲取1990、1998、2009年的山東省Landsat-5 TM影像、2018年的Landsat-8 OLI影像作為數據源；通過GEE平臺自動執(zhí)行影像拼接、輻射校正、去云等預處理。

(2)特征計算：所選特征包括纓帽變換亮度、綠度、濕度分量，歸一化植被指數(NDVI)，土壤調節(jié)植被指數(SAVI)，修正后的歸一化差異水體指數(MNDWI)，GLCM紋理特征。

圖2 基于GEE的農業(yè)大棚信息提取流程圖Fig.2 Flow chart of GEE-based agricultural greenhouse information extraction

(3)訓練集和測試集制作：首先計算各采樣點對應的特征值，并按照3∶1劃分訓練集和測試集。

(4)隨機森林分類：訓練隨機森林分類器，完成對驗證樣本的精度檢驗。

(5)制圖與分析：基于分類結果制作近30年山東省農業(yè)大棚空間分布專題圖，分析農業(yè)大棚時空動態(tài)變化的驅動力和影響因子。

2.2 特征計算

為了提高農業(yè)大棚的分類精度，本文基于纓帽變換等構建了一個多維特征空間，所有特征計算均在GEE平臺下編程實現。

2.2.1纓帽變換

纓帽變換(Tasseled-cap transformation，TCT)是將圖像的原始波段轉換為一組新的對植被較為敏感的波段。用于創(chuàng)建纓帽帶的系數是從圖像和經驗觀察中統(tǒng)計得出，其中第1個纓帽帶對應于圖像的整體亮度；第2個纓帽帶對應綠度；第3個纓帽帶通常被解釋為濕度(土壤或表面水分)。在本研究中，增加TCT特征可有效提高農業(yè)大棚與水澆地、露天植被之間的可分性。纓帽變換公式為

Y=cX+b

(1)

式中Y——變換后多光譜空間的像元矢量

X——變換前多光譜空間的像元矢量

c——變換矩陣b——常數

2.2.2光譜指數

除了使用TCT外，本文通過文獻調查和實驗比對選取歸一化植被指數(NDVI)、土壤調節(jié)植被指數(SAVI)、修正后的歸一化差異水體指數(MNDWI)3個光譜指數。其計算公式如下

(2)

(3)

(4)

式中ρ(G)——綠色波段反射率

ρ(R)——紅色波段反射率

ρ(N)——近紅外波段反射率

ρ(M)——短波紅外波段反射率

L——校正因子[21]

NDVI是用于提高植被和非植被的差異性，SAVI則能減小NDVI受土壤背景噪聲的影響。MNDWI指的是NDWI的修改版本[22]，由XU[23]提出。MNDWI是通過用MIR頻段代替NDWI的近紅外頻段得到的，具有更高的水體提取精度。

2.2.3紋理特征

根據實地采樣和目視判斷，發(fā)現山東省農業(yè)大棚具有明顯的幾何和紋理特征，如較寬的矩形邊界，外形一般有雙窗面、三折面、半拱圓形、連棟形等特定的外形特征。本研究選取6個最不相關的紋理特征，即均值(MEA)、標準偏差(STD)、同質性(HOM)、相異度(DIS)、熵(ENT)和角度秒矩(ASM)，用于區(qū)分大棚和其他地物之間差異。

選擇4個常用的移動窗口(3×3、5×5、7×7、11×11)對6個紋理特征進行比較，結果表明，各紋理特征均在3×3移動窗口中產生最高的分類精度。

2.3 分類類別

完成特征計算后，將Landsat中各波段的地表反射率、光譜和紋理特征添加到一個特征集中，以此作為特征數據集進行農業(yè)大棚分類器的訓練。土地利用類型分類類別如表1所示。

表1 分類類別Tab.1 Range of availability

本文采用半自動的采樣方法，首先通過實地調查、歷史樣本收集和目視解譯等方法，將得到的原始樣本在GEE云平臺上做聚類，將置信度、相似度高的作為樣本集，然后在GEE可視化界面中勾選樣本區(qū)域，所勾選的多邊形內的像素屬于同一類別，分別為農業(yè)大棚或其他土地利用類型。最后在勾選的1 300余個多邊形中隨機生成100 000個樣本點，每類土地利用類型的樣本點數量均為50 000個，確保覆蓋到典型的農業(yè)大棚和其他土地利用類型。其中隨機選取75 000個樣本點作為訓練集，剩余的25 000個作為測試集，用于測試混淆矩陣、總體分類精度等指標。

2.4 隨機森林算法

本研究使用隨機森林算法構建農業(yè)大棚分類模型。BREIMAN[24]在2001年提出了隨機森林算法，通過組合多個CART決策樹進行集成分類。在隨機森林訓練過程中，通過隨機選取2/3的訓練樣本和特征向量進行每棵決策樹的構建，模型輸出結果取決于所有決策樹的輸出，從而有效提高了模型的魯棒性和泛化能力，并擅長處理異常值和共線性變量。因此，隨機森林算法已被廣泛地用于遙感圖像分類領域。然而，隨機森林算法很少應用于農業(yè)大棚分類研究中，因此，本文從大空間尺度、長時間序列的角度驗證隨機森林算法對農業(yè)大棚的分類效果和泛化能力等。

隨機森林分類器的模型參數較少，比較關鍵的參數包括ntree和mtry。其中ntree是構成隨機森林的決策樹數量，mtry是每棵決策樹訓練時隨機抽取的特征數量。經過相關實驗，本文取ntree為150，以確保模型能夠收斂，mtry為特征總數的平方根。

2.5 精度評價

為了評價分類結果的質量，采用混淆矩陣和目視評估兩種方法。目視評估用于直接分析分類結果，是最直接、最基本的方法。混淆矩陣來自測試樣本，并計算總體分類精度、Kappa系數等指標。

3 結果與分析

3.1 農業(yè)大棚分類結果

圖3～5給出了2018年山東省濰坊市、臨沂市、棗莊市等局部典型地區(qū)的遙感影像和對應的農業(yè)大棚分類結果。可以看出，本文所提分類方法可以很好地提取農業(yè)大棚的空間分布信息。圖6給出了大棚分類結果的細節(jié)圖，可以看出分類結果與實際情況吻合程度較高，驗證了本文所提分類方法的性能。

圖3 濰坊市遙感影像及分類結果Fig.3 Remote sensing images and classification results of Weifang City

圖4 臨沂市遙感影像及分類結果Fig.4 Remote sensing images and classification results of Linyi City

圖5 棗莊市遙感影像及分類結果Fig.5 Remote sensing images and classification results of Zaozhuang City

圖6 農業(yè)大棚(鄉(xiāng)鎮(zhèn)級)遙感影像及分類結果Fig.6 Remote sensing images and classification results of agricultural greenhouses (township level)

同時，本文比較了不同機器學習方法在農業(yè)大棚提取中的性能，包括隨機森林算法(RF)、支持向量機(SVM)以及決策樹(DT)算法，使用相同的訓練樣本和驗證樣本以確保可比性和客觀性，所有實驗都是在GEE云平臺上完成。為了定量評估本文所提方法的準確性，表2給出了精度對比。由表2可知，RF性能優(yōu)于SVM和DT。RF的總體精度(OA)高于SVM，增幅為5.20～6.52個百分點。一個可能原因是在處理大規(guī)模、高維度數據集時RF比SVM效果更好。同樣RF的OA高于DT，增幅為6.09～8.44個百分點。這主要因為RF是多個決策樹的集合，與單個決策樹相比，RF具有更高的魯棒性，尤其可有效防止過擬合，有更高精度。

3.2 時空動態(tài)分析

本文得到的近30年山東省農業(yè)大棚的空間分布圖如圖7所示。從圖7可知，1990年農業(yè)大棚主要集中分布在壽光市及周邊。隨著山東省大力推廣農業(yè)大棚，全省各地如單縣、滕州市、蒼山縣、莘縣、萊西市等，借鑒壽光市農業(yè)大棚的成功案例，逐步引入了農業(yè)大棚。到2009年，以上市縣已經發(fā)展成為農業(yè)大棚區(qū)域的中心，并逐漸輻射到周邊市縣，初步形成多個具有較大規(guī)模的農業(yè)大棚聚集區(qū)。到2018年，全省除少部分受環(huán)境因素制約的地區(qū)外，其余市縣均在不同程度上實現了農業(yè)大棚的規(guī)模化建設，特別在魯中、魯西、魯西南等地區(qū)(如濰坊市、濟寧市、棗莊市、臨沂市、聊城市等地級市), 其中壽光市、昌樂縣、莘縣、萊西市等市縣農業(yè)大棚分布較集中，形成了多個規(guī)模龐大、功能齊全、布局完善的農業(yè)大棚示范區(qū)。

表2 RF、SVM和DT的精度比較Tab.2 Accuracy comparison of RF, SVM and DT

圖8 1990—2018年山東省農業(yè)大棚面積總量及面積增量、增速統(tǒng)計圖Fig.8 Diagrams of total area, area increment and growth rate of agricultural greenhouses in Shandong Province from 1990 to 2018

圖8為1990—2018年山東省農業(yè)大棚面積總量及增量、增速圖。由圖8可以得出，山東省農業(yè)大棚的發(fā)展經歷了緩慢-高速-緩慢的增長過程。1990年，全省農業(yè)大棚面積約為6.67 km2，1998年與1990年相比面積增量約為954.62 km2，面積平均增速為119.33 km2/a，面積達到961.29 km2。2009年全省農業(yè)大棚與1998年相比面積增量約為7 522.39 km2，平均增速達683.85 km2/a，總面積達8 483.68 km2。經歷了快速發(fā)展時期后，山東省積極響應中央深化供給側結構性改革，農業(yè)大棚從更快更廣向更精更專發(fā)展。到了2018年，與2009年相比增量約達1 435.72 km2，平均面積增速為159.52 km2/a，總面積達9 919.40 km2。經過30年的快速發(fā)展，山東省農業(yè)大棚從1990年的6.67 km2增加到2018年9 919.40 km2，平均增速為354.03 km2/a。

圖9為山東省1990—2018年農業(yè)大棚動態(tài)變化圖。從圖9可以看出，山東省農業(yè)大棚經過30年的快速發(fā)展，由最初空間分散、數量少，發(fā)展為規(guī)模龐大、功能完善的連片大棚聚集區(qū)。整體發(fā)展穩(wěn)中向好，面積穩(wěn)定增加。其中壽光市、滕州市為典型區(qū)域。

圖9 1990—2018年山東省農業(yè)大棚動態(tài)變化圖Fig.9 Dynamic changing map of agricultural greenhouses in Shandong Province from 1990 to 2018

20世紀80年代末，日光溫室(冬暖大棚)及黃瓜越冬栽培技術由遼南地區(qū)傳入山東地區(qū)[25]。1990年，壽光市較早發(fā)展農業(yè)大棚，總面積約為2.34 km2，經過30年的快速發(fā)展，始終處于全省農業(yè)大棚發(fā)展前列。到2018年，全市農業(yè)大棚約有862.92 km2。壽光市農業(yè)大棚起步早、基礎好，從1990年到2018年呈現數量穩(wěn)定增長、分布逐步外延、規(guī)模明顯增加的趨勢。不僅如此，壽光市更是積極帶動周邊市縣如青州市、昌樂縣、廣饒縣等快速發(fā)展，形成了全省面積最大、分布最廣、功能最完善的農業(yè)大棚聚集區(qū)。

全省其他市縣的農業(yè)大棚以滕州市為例，1990—1998年空間零星分布，數量少，規(guī)模小，主要分布于村莊周邊。1998—2009年得到快速發(fā)展，分布變廣且呈現進一步擴展趨勢，分布格局開始出現區(qū)域連片現象，發(fā)展初具規(guī)模，在3個時期中相對擴展比例最大，相對擴張速度最快，數量增長最多。2009—2018年，全省農業(yè)大棚基本空間分布格局已經形成，增長速度降緩，規(guī)模穩(wěn)步擴張。

圖10 1990—2018年山東省農業(yè)大棚2 km格網面積比例圖Fig.10 Proportion of 2 km grid area of agricultural greenhouses in Shandong Province from 1990 to 2018

為了進一步精準獲取農業(yè)大棚空間分布及動態(tài)變化，本文采用基于格網統(tǒng)計的方法將時空數據轉換到一個高分辨率基礎地理單元。本文通過實驗對比1、2、5、10、20 km格網后得出，2 km格網農業(yè)大棚空間分布的可視化效果最好。因此本文按2 km格網統(tǒng)計該區(qū)域農業(yè)大棚面積占比，共劃分為9等(圖10)，依次為1等(0～1%)、2等(1%～5%)、3等(5%～10%)、4等(10%～20%)、5等(20%～30%)、6等(30%～40%)、7等(40%～50%)、8等(50%～70%)、9等(70%以上)。如圖10所示，1990年全省農業(yè)大棚主要分布于壽光市，聚集度最高為5等(20%～30%)。隨著農業(yè)大棚的進一步推廣，在1998年，壽光市、青州市、昌樂縣出現了第1個9等農業(yè)大棚聚集區(qū)。進入21世紀以來，山東省大力扶持農業(yè)大棚產業(yè)，于2009年形成了6個9等農業(yè)大棚聚集區(qū)，例如以萊西市為中心，帶動平度市、即墨市等發(fā)展農業(yè)大棚；聊城市以東昌府區(qū)、莘縣為示范區(qū)，帶動陽谷縣、冠縣等周邊縣區(qū)農業(yè)大棚的發(fā)展；另外，滕州市、蒼山縣等利用自身發(fā)展優(yōu)勢，也積極向周邊市縣推廣農業(yè)大棚。到2018年，全省的農業(yè)大棚已經形成“一超多強”的發(fā)展態(tài)勢，共形成7個9等農業(yè)大棚聚集區(qū)，其中壽光市領跑全省，并進一步擴張。

通過以上分析可以得出，山東省農業(yè)大棚發(fā)展態(tài)勢良好，說明山東省在傳統(tǒng)農業(yè)向高效集約農業(yè)轉型卓有成效，為山東省由農業(yè)大省向農業(yè)強省的轉變奠定了堅實的基礎。

3.3 驅動力分析

分析農業(yè)大棚背后的驅動力和影響因子可以很好地總結發(fā)展經驗和制定未來規(guī)劃。通常，驅動力分析主要是從相互關聯的自然和社會因素角度進行研究[26-27]。在區(qū)域范圍內，社會因素對農業(yè)大棚的影響更為重要：

(1)國家政策因素：隨著“菜籃子工程”的實施，國家制定了詳細的幫扶政策來鼓勵農業(yè)大棚的建設。壽光三元朱村溫室的成功引起了各級政府的高度重視。1997年山東省農機局基于壽光成功案例部署了全省農業(yè)大棚工作[28]。從圖10可以看出，1990年到1998年期間，壽光市農業(yè)大棚得到快速發(fā)展，成為全省設施農業(yè)典型示范區(qū)。進入21世紀，山東省繼續(xù)完善農業(yè)大棚相關政策。在1998—2009年期間，各地積極學習壽光大棚經驗和技術，在全省范圍內出現了不同規(guī)模等級的大棚連片區(qū)，整體發(fā)展欣欣向榮。2013—2019年的中央一號文件都強調了發(fā)展現代設施農業(yè)的重要性。從圖10中可以看出，歷經30年發(fā)展，山東省農業(yè)大棚保持良好的發(fā)展態(tài)勢，形成了多個區(qū)域連片、輻射帶動效應明顯的農業(yè)大棚示范區(qū)。

(2)經濟與市場因素：人們對新鮮瓜果蔬菜的需求日益增長，尤其是在冬季，市場得不到及時供應，由此刺激了供給側的快速發(fā)展。從陳香等[29]研究中分析2000—2009年間山東省蔬菜進出口值，并結合圖10的動態(tài)變化可以看出，由于國內國際市場的需求量增大，不僅促進了山東省出現規(guī)模更大、功能更全的大棚代替中小棚，同時呈現出連片和聚集現象。隨著經濟全球化影響的進一步加深，山東省與韓國、日本等東北亞國家合作更加密切，發(fā)展設施蔬菜出口，可以促進生產和銷售的多元化。面對更加廣闊的海外市場，也給山東省的設施蔬菜提出了更高的要求和期望，由此進一步促進了山東省農業(yè)大棚的發(fā)展和優(yōu)化升級。

(3)生產技術因素：山東省農業(yè)大棚發(fā)展初期主要依靠人工操作，以經驗判斷為主，生產效率較低。隨著科技的高速發(fā)展，山東省農業(yè)大棚的生產水平日益提高，除了引入自動化、智能化控制系統(tǒng)外，各級政府還鼓勵科研及教育機構定向指導生產，發(fā)展關鍵共性技術，高效科學地指導生產。

4 結論

(1)通過利用GEE云平臺和多時相Landsat影像構建時序光譜特征和紋理特征，基于隨機森林分類器實現了山東省近30年農業(yè)大棚分類，并進行了時空動態(tài)分析。結果表明，本文方法取得了較高的分類精度，其平均總體精度為91.63%，實現了省域尺度的農業(yè)大棚空間分布提取。

(2)目前山東省農業(yè)大棚重點分布在魯中、魯西和魯西南等地區(qū)，如濰坊、濟寧、棗莊、臨沂、聊城等地級市, 其中壽光、昌樂、莘縣、平度等市縣分布集中。時空動態(tài)分析表明，山東省農業(yè)大棚在近30年間得到了快速發(fā)展。農業(yè)大棚面積從1990年的6.67 km2增加到2018年的9 919.40 km2，增長速度為354.03 km2/a。此外，1998年至2009年農業(yè)大棚增長速度最高，為683.85 km2/a，而1990年至1998年和2009年至2018年相對變化較小。