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

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

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

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

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

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

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

1 材料和方法

1.1 研究區(qū)概況

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

圖1 關中地區(qū)地理位置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 統(tǒng)計數據 冬小麥播種面積統(tǒng)計數據來源于陜西省統(tǒng)計年鑒[3]。本文收集的市級統(tǒng)計數據包括寶雞、咸陽、西安、銅川、渭南共5個市，時間跨度為2011—2017年?？h級統(tǒng)計數據包括寶雞8個市(縣、區(qū))，咸陽12個市(縣、區(qū))，西安6個市(縣、區(qū))，銅川3個市(縣、區(qū))，渭南11個市(縣、區(qū))，共計40個市(縣、區(qū))，時間跨度為2011—2016年。

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

1.3 研究方法

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

(1)

式中，a為常數?；贕EE平臺，以單個像元為單位，在每個像元位置上提取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時間變化規(guī)律來反映空間變化規(guī)律。提取從10月24日到12月19日的MODIS NDVI影像，建立NDVI值與時間的一元線性回歸模型，得到一幅含有兩個波段的影像，一個為斜率，另一個為截距。同理可得到次年4月15日到6月18日的一元線性回歸影像。

1.3.3 精度評價方法 在市級尺度和縣級尺度上，分別對冬小麥提取面積與農業(yè)統(tǒng)計面積進行相關性分析，評價提取精度時采用了決定系數(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與農業(yè)統(tǒng)計面積S′進行對比，面積提取精度P的計算公式為：

(4)

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

2 結果與分析

2.1 不同地物的NDVI時序特征

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 冬小麥種植面積時空變化規(guī)律

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

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

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

3 討論與結論

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

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

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

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