基于Google Earth Engine平臺的關中冬小麥面積時空變化監測

郭 新，王乃江，張玲玲，郭永強，褚曉升，馮 浩2,

(1.西北農林科技大學水利與建筑工程學院，陜西 楊凌 712100；2.西北農林科技大學水土保持研究所, 陜西 楊凌 712100；3.中國科學院水利部水土保持研究所，陜西 楊凌 712100)

小麥作為中國的主要糧食作物之一，其播種面積占全國糧食作物播種總面積的20%左右[1]。陜西省是我國重要的冬小麥生產基地[2]，關中地區冬小麥種植面積占陜西省冬小麥種植面積的83%[3]。因此，盡可能及時、準確地獲取關中地區冬小麥種植面積及其時空變化規律，對于發揮冬小麥的生產潛力，促進農業生產的可持續和高效發展，保障國家糧食安全具有至關重要的現實意義[4]。然而，傳統的冬小麥種植信息獲取方法耗時耗力，需要大量的經濟支出[5]，且傳統調查方法更新速度慢，易受人為因素影響[6]。

遙感技術因其時效性高、成本低、客觀性強、覆蓋范圍廣和信息量大等特點[7-8]，為大面積農作物信息提取和監測提供了有效的方法和手段[9]。由遙感數據反演而來的歸一化植被指數(normalized difference vegetation index，NDVI)能夠充分反映植被的季節性變化特征[10]，因而被廣泛應用于植被的生長動態監測、作物識別與土地覆蓋分類等研究領域。有學者利用冬小麥的物候特征，通過對關鍵生育期NDVI值或多個生育期NDVI變化值設定適當的閾值，建立冬小麥提取規則，進而獲取冬小麥的空間分布信息[11]。如楊小喚等[12]、李紅梅等[13]通過該方法分別提取了北京市和陜西關中地區的冬小麥種植面積。然而以上研究均是對幾個關鍵時相影像應用NDVI閾值法，卻忽略了冬小麥生育期內NDVI時序曲線的整體波動，且未考慮曲線的“峰”、“谷”特征。

以谷歌地球引擎(Google Earth Engine，GEE)為代表的云計算平臺，其公共數據庫提供了超過40 a的遙感影像[14]，且這些數據都經過預處理[15]。該平臺具有強大的數據處理能力，可進行復雜的空間可視化分析，目前國內外已經有基于該平臺的亞洲水稻種植面積提取研究[16]、全球城市土地監測[17]、土地利用變化研究[18-19]、新疆積雪覆蓋研究[20]、贛南柑橘果園分布制圖[21]等。該平臺為地物信息大范圍長時間序列的空間可視化分析提供了便利。

本文基于GEE平臺，根據關中地區冬小麥的物候特征，結合該地區典型地物NDVI時序曲線，采用NDVI重構增幅算法和光譜突變斜率，構建了冬小麥的提取模型，利用該模型識別了研究區內冬小麥種植面積，然后利用統計數據和地面調查數據分別從市級和縣級尺度對提取結果進行驗證，最終揭示關中地區冬小麥種植面積時空變化規律，以期為長時間序列大面積區域快速、準確地提取冬小麥種植面積提供一種普適性的模型和方法。

1 材料和方法

1.1 研究區概況

關中地區位于陜西省中部(圖1)，介于106°18′～110°38′E、33°35′～35°52′N之間，區域總面積5.56×104km2，包括寶雞、咸陽、西安、銅川、渭南5個省轄地級市。關中地區西起寶雞峽，東至黃河，北部為黃土高原南緣，中部分布著渭河沖積平原，南部是秦嶺北坡，因此該地區的地勢特點為南北高、中間低。海拔高度322～3 689 m。該地區屬于暖溫帶半濕潤大陸性季風氣候區，四季分明，多年平均氣溫7. 8～13.5℃，多年平均降水量500～800 mm。冬小麥是該地區主要的夏收作物，其生育期為10月至次年6月，常年種植面積約為9×105hm2。

圖1 關中地區地理位置Fig.1 Location of Guanzhong Region

1.2 數據來源及處理

1.2.1 MODIS NDVI數據 本文采用2010年10月—2017年6月冬小麥生育期內的Terra/MOD13Q1NDVI數據，中分辨率成像光譜儀(moderate resolution imaging spectroradiometer，MODIS)數據合成周期為16 d，空間分辨率為250 m，來源于GEE的公開數據庫。MODISNDVI數據有效值為[-2000,10000]，比例因子為0.0001。本文將NDVI值域轉化為[-0.2,1]，便于分析及確定分類閾值。本文對遙感影像數據的處理在GEE中完成。

1.2.2 統計數據 冬小麥播種面積統計數據來源于陜西省統計年鑒[3]。本文收集的市級統計數據包括寶雞、咸陽、西安、銅川、渭南共5個市，時間跨度為2011—2017年。縣級統計數據包括寶雞8個市(縣、區)，咸陽12個市(縣、區)，西安6個市(縣、區)，銅川3個市(縣、區)，渭南11個市(縣、區)，共計40個市(縣、區)，時間跨度為2011—2016年。

1.2.3 地面調查數據 地面調查數據通過實地采樣調查和目視解譯獲取。實地采樣調查時間為2017年12月9日和2019年3月15日，用手持式GPS定位儀獲取冬小麥樣點的經緯度信息。目視解譯則是利用Google Earth Pro的高分辨率影像對冬小麥種植區域的物候和紋理進行判別。地面調查樣點共有757個(圖1)，調查樣點均勻分布在關中地區各縣，每個調查樣點均選取常年種植且面積較大的冬小麥田塊，以保證純凈的冬小麥像元。地面調查數據隨機按7∶3比例分為2組，分別用于分類過程中訓練樣本的選取和分類結果的精度驗證。

1.3 研究方法

1.3.1NDVI重構增幅算法 冬小麥從播種到收獲，其葉片覆蓋度從“無”到“有”再到“無”的過程在NDVI曲線上表現出明顯的“峰”、“谷”形狀特征。本文采用NDVI重構增幅算法[22]，即計算冬小麥生育期內第一個“谷”到最高 “峰”的NDVI增長幅度(NDVIincrease)，其計算公式為：

(1)

式中，a為常數。基于GEE平臺，以單個像元為單位，在每個像元位置上提取9月15日至11月15日的NDVI最小值，得到一幅最小NDVI圖像，將其記為NDVImin；同理，在每個像元位置上提取次年3月1日至4月31日的NDVI最大值，得到一幅最大NDVI圖像，將其記為NDVImax。

1.3.2 光譜突變斜率 冬小麥的NDVI曲線隨時間呈現周期性的上升或下降，本文提出光譜突變斜率(Slope)這一概念，表示NDVI值隨時間的變化特征，其計算公式為：

(2)

(3)

式中,物理量右下角數字表示日期，NDVI和Date表示該日期的NDVI值和儒略日。b,c分別為常數。

基于GEE平臺，在每個像元位置上建立NDVI和時間的一元線性回歸模型，是以單個像元的NDVI時間變化規律來反映空間變化規律。提取從10月24日到12月19日的MODIS NDVI影像，建立NDVI值與時間的一元線性回歸模型，得到一幅含有兩個波段的影像，一個為斜率，另一個為截距。同理可得到次年4月15日到6月18日的一元線性回歸影像。

1.3.3 精度評價方法 在市級尺度和縣級尺度上，分別對冬小麥提取面積與農業統計面積進行相關性分析，評價提取精度時采用了決定系數(coefficient of determination，R2)均方根誤差(root mean square error，RMSE)和一致性指標(index of agreement，d)3個指標。其中R2越大，RMSE越小，表明提取精度越高；d≥0.9、0.8≤d<0.9、0.7≤d<0.8和d<0.7分別表示提取精度極高、高、一般和差[23-24]。

同時，為了定量分析本研究提取結果的精度，將提取的冬小麥面積S與農業統計面積S′進行對比，面積提取精度P的計算公式為：

(4)

式中，P值介于0～1之間，越接近1提取精度越高。

2 結果與分析

2.1 不同地物的NDVI時序特征

對研究區域內林地、建設用地、水體、果園和冬小麥等典型地物分別選取10個樣本區域，獲取樣本區域2011～2017的NDVI時序曲線，剔除異常曲線，各地物類型分別取平均值，再進行年際間的平均得到關中地區典型地物NDVI時序曲線(圖2)。從圖中可以看出，不同地物NDVI變化不盡相同，各有特點。林地的NDVI曲線在5—9月水平較高，波峰持續時間長，并且變化非常平緩。由于果樹的種植間距較大，覆蓋度不均勻，相對較為稀疏，因此果園的NDVI水平低于林地，但變化趨勢與林地基本一致。建設用地和水體的NDVI曲線整體水平很低且曲線波動較小，這是由于這些地物在近紅外波段上反射率較低，特別是水體，其NDVI值甚至會出現負值。另外，由于建設用地通常與城市綠地混合，所以建設用地的NDVI水平相對高于水體，且其趨勢與林地一致。

圖2 關中地區典型地物NDVI時序曲線圖Fig.2 NDVI time-series curves of typical objects in Guanzhong Region

冬小麥的NDVI曲線具有非常明顯的季節性變化特點。關中地區冬小麥一般在10月上旬播種，直到12月均處于緩慢生長階段，NDVI曲線呈上升趨勢，而此時其他植被基本進入越冬期，NDVI曲線呈下降趨勢。隨著氣溫降低，1—2月冬小麥基本停止生長。在2—4月冬小麥拔節～孕穗～抽穗期，冬小麥生長快速，NDVI迅速上升并達到峰值。冬小麥NDVI值從播種時的谷值到峰值有較大的增長幅度。之后隨著冬小麥逐漸灌漿成熟，NDVI曲線呈明顯下降趨勢，而此時其他植被正處于生長期，NDVI不斷增加。根據以上特點提取關中地區冬小麥種植面積。

2.2 冬小麥種植面積提取結果及其精度驗證

結合NDVI時序曲線特征、NDVI重構增幅算法和光譜突變斜率，由訓練樣本得知，a、b、c分別為0.49、0、-2.6，代入公式(1)～(3)中，提取2011—2017年關中地區冬小麥種植面積，提取結果見圖3。

由圖3可以看出，2011—2017年關中地區冬小麥空間分布雖然在年際間有所變化，但變化并不大。受地形影響，冬小麥種植區域主要集中分布在中部地勢平坦的關中平原，在北部及南部地勢較高的地區分布較少。冬小麥在寶雞、咸陽、西安和渭南分布較廣，是由于這4個市土壤肥沃、灌溉條件便利、自然條件較為優越,而在銅川市分布相對稀疏。從縣域行政區劃單元看，關中地區的冬小麥種植較大的區、縣主要分布在寶雞市的陳倉區、鳳翔縣、岐山縣、扶風縣，咸陽市的永壽縣、乾縣、武功縣、興平縣和涇陽縣，西安市的戶縣、藍田縣，以及渭南市的富平縣、蒲城縣、臨渭區、大荔縣等。以上地區形成了關中地區冬小麥空間集聚區。

圖3 2011-2017年關中地區冬小麥種植分布圖Fig.3 Distribution of winter wheat planting in Guanzhong Region from 2011 to 2017

分別以市和縣為統計單位，在研究區內對遙感提取冬小麥面積與統計數據中冬小麥播種面積進行相關性分析，檢驗冬小麥提取結果的精度。

圖4分析結果顯示，在市級和縣級水平，遙感提取面積和農業統計面積均呈現出極顯著的線性相關關系，R2分別為0.82和0.62，RMSE分別為42.98×103hm2和14.77×103hm2，d分別為0.95和0.84，市級和縣級的提取精度分別為“極高”、“高”，市級尺度擬合效果優于縣級尺度。可能因為在縣級尺度內正負誤差相互抵消，從而表現為在市級尺度上精度有所提高。即遙感數據的空間分辨率對提取結果產生了尺度效應[25]。

圖4 2011—2017年冬小麥遙感提取面積在市級(A)和縣級(B)水平與農業統計面積的對比Fig.4 Comparison between extracted winter wheat areas and statistical areas at city level (A) and county level (B), respectively, from 2011 to 2017

匯總關中地區各市多年冬小麥提取面積，與農業統計面積進行對比驗證，驗證結果的統計值見表1。可以看出，各市累計7 a中遙感提取面積與農業統計面積的一致性小于60%的概率僅為11.4%。其中渭南、咸陽的提取精度最高，7 a間面積一致性>80%年份達到5 a以上，且沒有出現<60%的年份；寶雞、西安的提取精度較高，>80%的年份均為4 a，且全部年份提取精度均在60%以上；銅川市的提取結果精度較低，僅有42.8%的年份提取精度>60%。這可能與銅川市地形地貌條件復雜有關，銅川地處黃土丘陵溝壑區，冬小麥種植區域分布較為破碎。當利用MODIS 250 m遙感影像識別冬小麥時，種植區域較大且集中的地塊可以比較容易被提取出來，而地塊面積小且分散破碎的地區則提取精度有限。另外，與關中地區其他市相比，銅川市冬小麥種植面積較小，故該市提取精度不高對整個地區影響不大。

表1 2011—2017年關中地區各市冬小麥種植 面積一致性累計值

另外，在研究區內隨機選取227個采樣點，驗證2017年遙感提取結果與地面調查數據的空間一致性。結果表明，227個冬小麥驗證樣點中，212個被正確分類，空間一致性精度為93.4%。

2.3 冬小麥種植面積時空變化規律

從市域行政區劃單元看，2011―2017年關中各市冬小麥種植面積變化趨勢不盡相同(圖5)。西安和寶雞的冬小麥種植面積變化趨勢較為一致，均呈下降趨勢。尤其是西安，整體收縮量達107×103hm2；渭南、咸陽和銅川的冬小麥播種面積有增有減，總體呈波動上升趨勢。2011―2017年扶風縣、咸陽市區、周至縣、戶縣、西安市區、長安區及渭南市區冬小麥種植面積呈明顯下降趨勢；永壽縣、乾縣和蒲城縣冬小麥種植面積明顯增加。

就關中地區整體而言，2011―2017年冬小麥種植面積呈下降趨勢，由2011年的981.96×103hm2縮減至2017年的898.74×103hm2，減少了83.22×103hm2(8.47%)。從空間分布上看，關中地區2011―2017年冬小麥種植面積呈“北增南減”的時空變化格局：其中新增冬小麥種植面積221.40×103hm2，分布在永壽縣、乾縣和蒲城縣，從Google Earth Pro上可以看到部分縣由林地向耕地的轉化，也證明各地土地整治工作卓有成效；其中縮減冬小麥種植面積304.62×103hm2，主要分布在扶風縣、西安市區、長安區、周至縣和戶縣、咸陽市區和渭南市區。

經過調查發現，關中地區冬小麥種植面積縮減的主要原因有以下3點：(1)農業收入低、投入高是該地區農戶棄耕的根本原因。其中，地下水超采導致的灌溉成本的提高及青壯年勞動力外出打工引發的勞動成本的提高可能是主導因素；(2)由于當地農業產業調整，農戶因地制宜發展苗圃、果樹等高收益作物；(3)隨著城市化進程的加快，犧牲大量的耕地轉化為建設用地是冬小麥種植面積減少的主要原因。

3 討論與結論

本文根據關中地區典型地物NDVI時序曲線變化特征，結合冬小麥物候特點，基于NDVI重構增幅算法和光譜突變斜率，且利用GEE強大的數據分析和數據處理能力，提取了2011―2017年的冬小麥種植區域，并檢驗了遙感提取結果的精度，揭示了2011―2017年關中地區冬小麥時空變化規律。該方法綜合考慮冬小麥NDVI時序曲線，參數少，能夠快速準確地提取冬小麥種植區域，有較高的實用價值和現實意義。主要結論如下：

(1)利用統計年鑒對冬小麥提取結果的驗證表明，在市級和縣級尺度R2分別為0.82和0.62，d分別為0.95和0.84；提取結果與實地調查數據的空間一致性精度為93.4%，因此本文所使用的提取方法可滿足大區域冬小麥面積提取的需求精度。

(2)關中地區冬小麥主要分布在中部較為平坦的關中平原，在北部及南部地勢較高的地區分布較少。

(3)2011―2017年，關中地區冬小麥種植面積呈“北增南減”的時空變化格局：乾縣、永壽縣和蒲城縣等地冬小麥種植面積增加；扶風縣、長安區、周至縣、戶縣、咸陽市區和渭南市區等地冬小麥種植面積減少。整體看來，關中地區冬小麥種植面積呈下降趨勢，減少了83.22×103hm2(8.47%)。銅川市提取精度較低，因為使用250 m空間分辨率的遙感影像導致地形復雜、冬小麥分布破碎的種植區不易被識別。因此，選擇分辨率更高的遙感影像、混合像元分解等問題將是后續研究的重點。