基于ANN多源遙感數據融合的耕地種植強度 估算方法研究

徐 猛，陶建斌，吳琪凡

（地理過程分析與模擬湖北省重點實驗室/華中師范大學城市與環境科學學院，湖北武漢430079）

0 引言

農業生產是人類社會存在和發展的基石，耕地是進行農業生產的基本資源和條件。隨著我國城市化迅速發展，每年約有20萬hm2耕地被轉化為建設用地［1-3］。同時，隨著人口持續增長、飲食結構改變和非食物用途（如生物能源）農產品消耗的增加，人類對農產品的需求不斷增長。2015年末全國人口為13.75億，預計到2030年全國人口將增加至14.5億［4］，增長5.45%，而同時期糧食需求將增長47.19%，達到9.17億t［5］。如何緩解日益加劇的人地矛盾，是國際社會普遍關注的焦點問題。在耕地面積減少、糧食需求持續增長和部分耕地質量惡化等多重壓力下［6］，耕地集約利用成為保障國家糧食安全的重要選項。同時，由于城市化快速發展對農業勞動力的虹吸效應，以及市場經濟的發展對農民種糧收益的沖擊，部分地區農作物種植頻率呈現下降趨勢（表現為耕地撂荒和冬閑田）。因此，耕地的集約利用成為解決未來世界糧食安全的重要途徑，引起學界和業界的廣泛關注。

當前，國內外研究者多數將耕地的集約利用定義為種植頻率（單季和雙季等）［7-13］或者復種指數［14-20］，也有部分研究者從投入產出、生產管理和技術角度來對耕地利用強度進行定義［21-27］。Kehoe等［23］從輸入、輸出、系統3個測度總結了現有的對Agricultural Land-use Intensity（簡稱ALI）的定義。Jiang等［24］認為ALI通常有3種測度：種植頻率、農業產出、農業投入。然而，無論是種植頻率還是復種指數，都是對種植模式的硬性劃分（如單季或雙季）。一方面，盡管各種復種指數定義不盡相同，但本質上都是種植頻率在空間尺度上的聚合或是時間尺度上的綜合。另一方面，基于統計數據得到的一系列耕地指標缺乏詳細的空間分布信息，忽略了行政區內部耕地集約利用水平的空間異質性。因此現有評價體系存在時、空分辨率上的矛盾，不能滿足耕地集約化利用的精細評價的要求。

在國外，很早就有了Cropping Intensity或ALI的概念［28-30］，但卻有許多不同的定義。在對Cropping Intensity進行系統梳理的基礎上，文章定義種植強度指數對耕地集約化利用進行精確表征，通過構建人工神經網絡（Artificial Neural Network，ANN）對時間序列MODIS數據與種植強度間的映射關系進行表達，以湖北省為研究區進行種植強度估算。該文研究方法可獲取高時空分辨率的種植強度數據集，對推動作物制圖的進一步發展，提高耕地集約化精細監測具有重要意義，為全球變化過程的分析與模擬提供理論、方法和數據支撐。

1 研究區及數據

1.1 研究區

湖北省地處我國中部（東經108°21′42″～116°07′50″、北緯29°01′53″～33°6′47″），面積18.59萬km2，約占全國國土面積的1.97%，耕地面積占到全國耕地面積的3.89%。湖北省東、西、北三面環山，中部低平。除部分山區外，大部分地區屬亞熱帶季風氣候，光熱充足，降水豐沛。優越的自然環境和悠久的種植傳統使湖北省成為著名的魚米之鄉，其糧油生產在全國占有重要地位。

1.2 數據來源及數據預處理

該文采用的主要數據包括：（1）MOD13Q1（MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid）250 m分辨率16 d合成的植被指數數據［31］；（2）Landsat 8遙感影像［32］。另外還采用了其他數據作輔助，包括Google Earth歷史影像，GlobeLand30 —2010全球30 m空間分辨率土地覆蓋數據［33］，中國1∶400萬矢量數據集。

MOD13Q1是陸地植被指數數據，具有250 m的空間分辨率和16 d的時間分辨率，在地表植被監測領域的應用十分廣泛。MOD13Q1數據包括增強型植被指數（Enhanced Vegetation Index，EVI）、歸一化植被指數（Normal Differential Vegetation Index，NDVI）、像素質量圖層（Pixel Reliability，PR）等12個波段。該文采用2015年MOD13Q1數據，研究區所涉及圖幅編號為h27v05、h27v06、h28v05、h28v06。高分辨率遙感數據采用了3期Landsat 8多光譜遙感數據（表1），空間分辨率為30 m，重返周期為16 d。該研究中，MOD13Q1數據用于神經網絡的輸入特征，Landsat 8遙感數據主要用于訓練樣本和驗證樣本的制作。

圖1 研究區地理區位及地物覆蓋類型（GlobeLand30—2010）Fig.1 The location of the study area and the land-cover types（GlobeLand30 —2010）

表1 Landsat 8實驗數據詳細信息Table 1 Landsat 8 experimental data details

數據預處理包括3個部分：（1）利用MODIS批量處理工具（MRT）對MOD13Q1數據進行預處理，重采樣至240 m，投影坐標選擇WGS_1984_UTM_Zone_49N，提取出EVI圖層，并以Landsat 8影像為基準對MOD13Q1數據進行配準。（2）MOD13Q1數據受到云、雪、陰影等諸多因素的影響，原始EVI數據存在很多噪聲，需要通過數學方法對原始EVI數據進行重建。TIMESAT是用于平滑長時間序列EVI數據的有效工具［34］，提供多種數據重建方法，該文采用Savitsky-Golay算法（簡稱S-G濾波算法），這是一種基于最小二乘法的移動窗口加權平均的數據重建方法［35］。S-G濾波算法在去除噪聲的同時，能夠清晰的保留局部細節，并且對EVI曲線的峰值和峰寬保真度較高［36］。（3）提取耕地掩膜，對GlobalLand 30 —2010年數據產品進行投影變換、裁剪、重分類等處理，提取湖北省的耕地掩膜。

2 研究方法

利用神經網絡估算耕地種植強度的基本假設是：多時相遙感影像的EVI數據反映了耕地的利用狀態，耕地的不同利用狀態與耕地種植強度之間存在映射關系。該文將MOD13Q1時間序列EVI數據和基于Landsat 8數據獲取的樣本區種植強度分別作為神經網絡的輸入和輸出層，訓練網絡模型，利用該網絡模型估算整體耕地種植強度。

2.1 種植強度指數

種植強度指數是指在亞像元尺度下綜合表征農作物的空間分布和種植頻率的指數，其公式為：

式（1）中，CII為種植強度指數，Dcrop為1個MODIS像元內Landsat 8遙感影像上雙季作物的像元個數，Scrop為1個MODIS像元內Landsat 8遙感影像上單季作物的像元個數，n為1個MODIS像元內Landsat 8遙感影像上的像元總數。

CII為MODIS像元內的種植頻率的均值，值域范圍為0～1（圖2）。在地塊完整的雙季作物種植區，其種植強度為1或接近于1。在地塊完整的單季作物種植區，其種植強度為0.5或接近0.5。非耕地區域的種植強度則為0（忽略研究區內少量的三季作物）。

圖2 耕地種植強度示意圖Fig.2 Diagram of cropping intensity

2.2 樣本種植強度

Landsat 8遙感影像用于樣本區耕地種植頻率的分類。在種植頻率分類結果的基礎上采用ArcGIS的分區統計與重采樣工具實現種植頻率的概化工作。

湖北省耕地種植模式主要以單季和雙季為主［37］。單雙季作物的物候期存在較大差異，單季作物生長起始期較晚，大概在每年4月初。雙季作物第1個生長季約在每年1月中下旬開始，在4月中下旬結束，EVI曲線在第129、145 d處明顯下降（圖3）。以2015年為例，在1—4月的Landsat 8真彩色合成圖像上，單季作物種植區呈現暗灰色，雙季作物種植區呈現暗綠色（圖4）。選取MOD13Q1的001、145這2期EVI影像，將001、145、001分別賦予紅、綠、藍通道，能夠明顯地凸顯出雙季作物種植區，根據加色法原理，雙季作物種植區呈現出品紅色（圖5）。借此結合Google Earth高清影像以輔助判讀Landsat 8遙感影像的非監督分類結果。樣本區的選擇應包括耕地、水體、天然植被、建設用地、裸地等多種地物類型，且面積不應小于研究區面積的15%，以增強模型的魯棒性和泛化性。該文采用非監督分類ISODATA分類方法對樣本區Landsat 8遙感影像進行分類，初始類數量為20類，迭代次數為15次。然后將初始分類結果重分類為雙季作物、單季作物和其它3類。借助ArcGIS創建240 m的格網，分區統計每個網格單元內種植頻率的平均值，得到每個網格單元的種植強度，值域為0～1（圖6）。

圖3 湖北省耕地單季作物與雙季作物時間序列EVI曲線Fig.3 EVI profiles of single cropping crops and double cropping crops of cropland in Hubei Province

圖4 2015年湖北省耕地單季作物與雙季作物典型分布區Fig.4 Typical distribution areas of single cropping crops and double cropping crops of cropland in Hubei Province in 2015

圖5 2015年湖北省耕地MOD13Q1 EVI假彩色合成Fig.5 False color composite of MOD13Q1 EVI of cropland in Hubei Province in 2015

圖6 基于Landsat 8影像的訓練樣本區數據：a. 種植頻率；b. 耕地種植強度Fig.6 Cropping frequency and intensity of training sample area based on Landsat 8

2.3 模型訓練

人工神經網絡是一種模擬生物神經網絡，進行分布式并行信息處理的數學模型，它具有自主學習的能力，通過不斷改變神經節點之間的連接權重，完成由輸入數據到目標數據的復雜映射［38］。BP神經網絡又稱為前饋型神經網絡，是人工神經網絡中應用最為廣泛的模型之一，其由3部分組成：輸入層、隱含層（可以是單層，也可以是多層）、輸出層。BP神經網絡通過誤差的反向傳播，不斷調整輸入層與隱含層、隱含層與輸出層之間的權重和誤差，使網絡輸出與實際輸出的誤差最小化［39］。該文利用MATLAB設計BP神經網絡自主學習樣本區23期EVI數據集與耕地種植強度間的映射關系，樣本區23期EVI作為輸入層，耕地種植強度作為輸出層。在隱含層數量及隱含層神經元個數的選擇上，該文參考Kolomogorov定理。Hecht-Nielsen指出，映射網絡存在Kolomogorov定理，即對于任意連續函數f ：［0，1］n→Rm都可以用1個3層前向神經網絡實現，且中間層神經元個數不多于2n+1個［40］。因此該文設計的BP神經網絡為3層：輸入層、輸出層和隱含層。關于神經元個數的選擇，在遵循Kolomogorov定理的基礎上，通過多次試驗，最終將神經元個數定為15，即構成1個23-15-1的BP神經網絡模型，見圖7。神經網絡其他參數設置如下：訓練函數TRAINLM，學習函數LEARNGDM，隱含層與輸出層傳遞函數均為LOGSIG。

圖7 BP神經網絡模型Fig.7 Back propagation neural network model

3 結果與分析

3.1 樣本精度驗證

基于神經網絡估算的耕地種植強度其準確性受樣本精度影響，為了利用驗證樣本對估算結果進行精度評價，需要保證樣本精度。該文從Google Earth歷史影像（2015年3月6日、2015年3月25日、2015年3月27日、2015年10月11日）上選取了其他、單季作物、雙季作物3類共計1 218個真實樣本（圖8），以驗證樣本精度。結果表明，基于Landsat 8影像得到的樣本總體精度達到了93.51%，Kappa系數為0.914（表2），證明了樣本具有很高的可信度。

圖8 2015年湖北省耕地真實樣本空間分布Fig.8 Spatial distribution of truth sample of cropland in Hubei Province in 2015

表2 樣本分類結果混淆矩陣Table 2 Confusion matrix of sample classification results

3.2 種植強度精度驗證

該文通過種植強度指數和人工神經網絡估算方法對驗證樣本區進行處理，制作驗證樣本區耕地種植強度，并與BP神經網絡提取的種植強度結果作對比。從目視結果來看，Landsat 8遙感影像上的種植頻率分布與BP神經網絡提取的耕地種植強度呈現出較高的一致性（圖9）。

利用回歸分析驗證BP神經網絡提取的湖北省耕地種植強度的精度。圖10顯示了驗證樣本區基于Landsat 8影像的耕地種植強度與BP神經網絡提取的耕地種植強度的散點圖和線性關系。R2達到0.923，回歸方程系數為0.948，截距項為0.023，說明二者之間的數據偏差很小，驗證了BP神經網絡提取耕地種植強度的可靠性。

圖9 a. 基于Landsat 8影像分類結果的驗證樣本區種植頻率 b. 基于BP神經網絡提取的驗證樣本區耕地種植強度Fig.9 a. Cropping frequency in validation sample area b. Cropping intensity in validation sample area extracted from BP neural network

圖10 驗證樣本區耕地種植強度與BP神經網絡耕地種植強度散點圖Fig.10 Scatter plot of cropping intensity in validation sample area extracted from Landsat 8 and cropping intensity in validation sample area extracted from BP neural network

3.3 湖北省耕地種植強度分析

該文利用神經網絡提取了2015年湖北省耕地種植強度（圖11）。基于BP神經網絡提取的種植強度實現了對耕地集約利用程度的軟劃分。BP神經網絡提取的種植強度與陶建斌等［41］提取的2015年湖北省種植頻率分布基本吻合（圖12）。圖13為陶建斌提取的2015年湖北省耕地種植頻率與該文BP神經網絡提取的耕地種植強度在江漢平原地區的對比，結果顯示BP神經網絡提取的耕地種植強度空間分辨率更高，細節表現能力更強，層次更加分明。

圖11 2015年湖北省耕地種植強度Fig.11 Cropping intensity in Hubei Province in 2015

圖12 陶建斌等提取的2015年湖北省耕地種植頻率［41］Fig.12 Cropping frequency from Tao Jianbin，et al in Hubei province in 2015

基于統計數據［42］計算得到2015年湖北省縣級尺度的耕地復種指數（圖14）。BP神經網絡提取的耕地種植強度與基于統計數據得到的耕地復種指數間總體上較為吻合，復種指數較高的江漢平原地區對應著種植強度高值區，復種指數較低的西部山區對應著種植強度的低值區。然而相較之下，BP神經網絡提取的耕地種植強度分辨率更高，對耕地集約利用水平的評價更為精確。

圖13 2015年湖北省耕地種植頻率與種植強度對比a. 該文提取的種植頻率 b. 陶建斌等提取的種植強度Fig.13 Comparison of cropping frequency and cropping intensity in Hubei province in 2015

圖14 2015年湖北省耕地復種指數Fig.14 Multi-cropping index in Hubei Province in 2015

4 結論與討論

該文提出了種植強度的概念，并定義了種植強度指數以精確表征耕地的集約化利用。構建BP神經網絡模擬時間序列EVI數據集與種植強度間的映射關系。估算結果精度達到92.3%，證明了BP神經網絡提取耕地種植強度的可靠性。該方法充分發揮遙感技術多尺度的時空表達能力和地物時間譜特征捕獲能力的優勢，并結合人工神經網絡自主學習能力的特點，可獲得高時空融合的耕地種植強度數據集。

傳統的基于統計數據的耕地集約化研究，盡管在時間尺度上能夠保證良好的連續性，但基于行政區單元的統計數據往往忽略了內部的空間異質性，空間分辨率較低；而國內外現有的基于遙感數據的耕地集約化利用研究雖然保證了較高的空間精度，但在方法上大都采用了連續年份的平均；該文提出的耕地種植強度彌補了現有指標的不足，將耕地集約利用研究尺度縮小至亞像元級別，具有融合高時間分辨率和高空間分辨率的特點，能夠實現對耕地種植強度的精細監測。

該文的研究重點在于方法的探究，未來將關注該方法應用于更大區域尺度耕地種植強度的時空變化以及種植強度估算中的地域分異問題。同時，該文在使用BP神經網絡時未探究其作用機理，導致網絡模型訓練的不確定性缺乏合理解釋。