潘 力，夏浩銘,2，王瑞萌，牛文輝，田海峰，秦耀辰,2

基于Google Earth Engine的淮河流域越冬作物種植面積制圖

潘 力1，夏浩銘1,2※，王瑞萌1，牛文輝1，田海峰1，秦耀辰1,2

（1. 河南大學地理與環(huán)境學院/河南省地球系統(tǒng)觀測與模擬重點實驗室/黃河文明與可持續(xù)發(fā)展研究中心暨黃河文明省部共建協(xié)同創(chuàng)新中心，開封 475001；2. 黃河中下游數(shù)字地理技術教育部重點實驗室（河南大學），開封 475001）

越冬作物是中國重要的作物類型，其面積的變化不僅對中國糧食產(chǎn)量和經(jīng)濟產(chǎn)生直接的影響，還潛在地影響中國的糧食安全，因此有必要準確繪制越冬作物種植面積圖來為決策制定者提供科學參考。該研究以淮河流域為例，基于Google Earth Engine云平臺，融合時間序列Landsat-7/8和Sentinel-2A/B衛(wèi)星影像，采用S-G（Savitzky-Golay）濾波器重構作物時間序列的歸一化差異植被指數(shù)（Normalized Difference Vegetation Index，NDVI），根據(jù)不同植被類型物候期的差異，選取越冬作物生長旺盛期NDVI最大值、越冬作物播種期和收獲期中相應的NDVI最小值和中位數(shù)，在像元尺度上構建越冬作物提取算法，繪制淮河流域越冬作物的種植面積。研究結(jié)果表明，所構建算法能夠精確提取淮河流域越冬作物的種植面積，總體精度為95.8%，Kappa系數(shù)為0.912，該研究可為作物面積的提取和監(jiān)測提供方法參考。

遙感；作物；制圖；種植面積；Google Earth Engine；淮河流域

0 引 言

越冬作物為秋季播種，幼苗經(jīng)過冬季，并在第二年春季或夏季收割的農(nóng)作物，例如冬小麥、冬大蒜和冬油菜等，是中國北部重要的糧食作物，及時準確地獲取越冬作物種植面積對于預測糧食產(chǎn)量、價格和國家的糧食安全至關重要[1-3]，并可用于提高農(nóng)業(yè)生產(chǎn)力和韌性[4]。精準的區(qū)域作物種植分布圖為作物生長監(jiān)測和產(chǎn)量預測提供了基礎數(shù)據(jù)，有助于幫助決策者及生產(chǎn)者制定合理的政策和風險管理策略[5-6]。

通常作物種植面積的獲取方法包括農(nóng)業(yè)調(diào)查法和遙感分類法。傳統(tǒng)的農(nóng)業(yè)調(diào)查法包括實地問卷調(diào)查和訪談等，逐層統(tǒng)計上報，此類方法既耗時又費力，缺少精確的作物空間分布信息，且時效性較差。遙感分類法能夠基于一幅或多幅影像精確地進行作物分類，是一種非常有效的植被監(jiān)測方法，已廣泛應用于區(qū)域、國家和全球作物制圖等[7-9]。隨著遙感數(shù)據(jù)的激增，免費共享和遙感大數(shù)據(jù)處理方法的發(fā)展，利用遙感技術進行作物播種面積制圖逐漸成為作物種植數(shù)據(jù)統(tǒng)計的重要途徑。

作物種植面積遙感制圖常用的方法是利用衛(wèi)星數(shù)據(jù)和分類算法在地表反射特征和作物生長特征之間建立關系[10]。由于中分辨率成像光譜儀（Moderate-resolution Imaging Spectroradiometer，MODIS）具有高時間分辨率的特點，且具有很好的時間序列特征，目前已開發(fā)了多種基于MODIS數(shù)據(jù)的算法來監(jiān)測中等空間分辨率下越冬作物的空間分布、估產(chǎn)、時空變化以及物候特征等，如周亮等[11]使用MODIS數(shù)據(jù)構建了基于卷積神經(jīng)網(wǎng)絡的冬小麥估產(chǎn)模型，張錦水等[12]基于MODIS數(shù)據(jù)構建了2017—2019年冬小麥時間序列數(shù)據(jù)對河北省冬小麥空間分布進行識別，Ren等[13]基于MODIS數(shù)據(jù)對黃淮海平原2001—2016年冬小麥進行估產(chǎn)，黃健熙等[14]利用MODIS的葉面積指數(shù)產(chǎn)品提取了冬小麥返青期、拔節(jié)期、抽穗期和開花期。盡管MODIS數(shù)據(jù)在進行作物動態(tài)監(jiān)測上得到了很好的應用，但因受到空間分辨率的限制，難以滿足破碎化、異質(zhì)性小塊農(nóng)田的監(jiān)測需求。隨著免費的中高空間分辨率衛(wèi)星的迅速發(fā)展，Landsat影像（30 m）逐漸應用于遙感制圖，如Phalke等[15]利用160 000幅Landsat影像繪制了30 m空間分辨率的歐洲、中東、俄羅斯和中亞大部分地區(qū)的農(nóng)田地圖，Wang等[16]利用Landsat影像繪制了1999—2018年美國中西部30 m空間分辨率的玉米和大豆地圖。Landsat影像解決了MODIS影像空間分辨率不足的問題，但由于Landsat衛(wèi)星的時間分辨率較低，且影像經(jīng)常受到云層的干擾，導致很難使用Landsat影像構建作物完整平滑的時間序列曲線。為了減少因Landsat低時間分辨率而導致的時間序列曲線值缺失的影響，許多研究開始結(jié)合使用多個傳感器的融合影像來進行研究[17-18]，如Li等[19]結(jié)合MODIS和Landsat數(shù)據(jù)繪制了安徽宣城作物種類和輪作的空間分布圖，Liu等[8]結(jié)合Landsat-7/8和Sentinel-2數(shù)據(jù)繪制了中國30 m空間分辨率的種植強度圖，Lenco等[20]結(jié)合Sentinel-1和Sentinel-2在10 m空間分辨率下對土地覆被進行分類。多個傳感器數(shù)據(jù)源的整合可以提高影像的時間分辨率，減少云覆蓋等因素的干擾，進而構建作物完整生命周期曲線，利用不同時期的光譜差異來區(qū)分作物的類型。

綜上，本研究基于Google Earth Engine云平臺，通過融合Landsat-7/8和Sentinel-2A/B數(shù)據(jù)集，首先構建了30 m高時空分辨率時間序列影像數(shù)據(jù)集，并計算歸一化植被指數(shù)（Normalized Difference Vegetation Index，NDVI）；其次采用S-G（Savitzky-Golay）濾波器對NDVI時間序列進行平滑，基于平滑后的NDVI時間序列確定越冬作物生長旺盛期、越冬作物播種期和收獲期，并分別計算生長旺盛期NDVI最大值、播種期和收獲期相應的NDVI最小值和中位值；最后通過構建分類模型，在像元尺度上繪制2018年淮河流域越冬作物種植面積圖。本研究可為作物面積的提取和監(jiān)測提供科學依據(jù)和方法參考，為作物管理和產(chǎn)量預測提供數(shù)據(jù)支撐。

1 材料與方法

1.1 研究區(qū)概況

淮河流域位于中國東部南北氣候過渡帶（111°55′E～121°20′E，30°55′N～36°20′N），橫跨河南、安徽、江蘇、山東、湖北5 ?。▓D1）?；春恿饔蛭鞑?、西南部及東北部為山區(qū)和丘陵區(qū)，其余為平原。流域面積為27萬km2，山區(qū)約占總面積的1/3，平原約占總面積的2/3?；春恿饔蚴侵袊匾霓r(nóng)業(yè)種植區(qū)域，廣闊的平原為作物提供了適宜的生長環(huán)境，耕地面積占整個流域面積的68%，占全國耕地面積的12%，糧食產(chǎn)量約占全國總產(chǎn)量的17%[21]?；春恿饔蚍N植模式多為單季種植和兩季種植，其復種比例達到58.4%[22]。單季種植耕地多為單季水稻、單季玉米等，兩季種植耕地多為冬小麥-玉米、冬小麥-花生等。近年來，淮河流域耕地面積呈遞減趨勢，越冬作物作為淮河流域主要作物種植類型，為保障國家糧食產(chǎn)量和糧食安全，探明淮河流域現(xiàn)有越冬作物種植狀況對越冬作物未來生產(chǎn)潛力的評估和提升十分重要。

1.2 數(shù)據(jù)來源

1.2.1 Landsat-7/8和Sentinel-2衛(wèi)星數(shù)據(jù)及處理

Landsat-7/8衛(wèi)星攜帶增強型專題制圖儀（Enhanced Thematic Mapper，ETM+）和陸地成像儀（Operational Land Image ，OLI），空間分辨率為30 m，重訪周期為15 d[23]。Landsat-7 ETM+和Landsat-8 OLI影像分別對應“LANDSAT / LE07 / C01 / T1_TOA”和“LANDSAT / LC08 / C01 / T1_TOA”。Sentinel-2由2顆衛(wèi)星組成，攜帶多光譜成像儀（Multi-Spectral Instrument，MSI)，空間分辨率為10、20和60 m，單顆星重訪周期為10 d[24]。Sentinel-2 A/B對應“COPERNICUS / S2”。本研究以越冬作物的一個完整物候期為例，選取時間范圍為2017年9月1日至2018年6月30日，合計4 982景影像。

基于Google Earth Engine云平臺，利用JavaScript編程獲取2017年9月1日至2018年6月30日淮河流域逐像元Landsat-7/8和Sentinel-2A/B衛(wèi)星影像的總觀測影像幅數(shù)如圖2a和圖2b所示。同時，采用Google Earth Engine云平臺的CFMask算法來識別云、云陰影、卷云和冰/雪覆蓋，去除所有質(zhì)量差的觀察結(jié)果[25]，獲得研究區(qū)逐像元高質(zhì)量觀測影像幅數(shù)如圖2c和圖2d所示。結(jié)果表明，在研究期間，研究區(qū)單個像元上觀測影像幅數(shù)最少為60幅，最高為534幅，其中高質(zhì)量觀測影像幅數(shù)最少為12幅，最高為311幅。

1.2.2 樣本數(shù)據(jù)

為確保足夠數(shù)量和質(zhì)量的樣本點數(shù)據(jù)用于分類及精度驗證，本研究基于中國氣象網(wǎng)（http://data.cma.cn/）提供的觀測數(shù)據(jù)資料、Google Earth高分辨率影像和野外調(diào)查數(shù)據(jù)（2018年3月15日—2018年7月1日）獲取了共621個具有代表性的、典型的像元作為樣本點（圖 1），分別為越冬作物樣本點250個、春播作物樣本點87個、落葉森林樣本點76個、常綠森林樣本點67個、建筑用地樣本點120個和水體樣本點21個。本研究按照3∶1的比例將樣本劃分為訓練樣本（463個）和驗證樣本（158個）。

Google Earth影像（空間分辨率為1 m）能清晰地顯示耕地、河流、道路、房屋村莊等地表特征，可作為本研究制圖結(jié)果的直接驗證數(shù)據(jù)。同時，本研究還使用2019年河南省統(tǒng)計年鑒[26]來間接驗證本研究的算法。由于淮河流域橫跨河南省、山東省、安徽省、江蘇省和湖北省5省，每個省統(tǒng)計年鑒的統(tǒng)計指標并不一致，僅有河南省統(tǒng)計年鑒細化至縣市級別，因此本研究選取淮河流域所覆蓋的河南省62個縣市的統(tǒng)計數(shù)據(jù)用于驗證，由于其中9個縣市位于淮河流域的邊界，不完全屬于淮河流域，故剔除這9個縣市，最后保留河南省53個縣市作為驗證區(qū)域。

1.2.3 坡度數(shù)據(jù)

本研究選取數(shù)字高程模型（Digital Elevation Model，DEM）數(shù)據(jù)，計算研究區(qū)的坡度（°），DEM數(shù)據(jù)可在Google Earth Engine云平臺中免費獲取，空間分辨率為30 m。

1.3 研究方法

本研究的主要思路如圖3所示，主要包括4個方面：1）數(shù)據(jù)集構建，包括構建30 m影像數(shù)據(jù)集、時間序列NDVI計算、插值和重構；2）基于平滑后的時間序列NDVI數(shù)據(jù)集，結(jié)合越冬作物物候期構建越冬作物提取算法；3）逐像元繪制2018年淮河流域越冬作物分布圖；4）利用Google Earth高分辨率影像和統(tǒng)計數(shù)據(jù)驗證越冬作物提取算法精度。

1.3.1 NDVI時間序列數(shù)據(jù)集構建及處理

基于NDVI時間序列曲線識別作物是常見且有效的方法，NDVI的計算如（1）所示。

由于不同傳感器所獲得影像的時間分辨率不同，掃描區(qū)域存在重疊，為了獲取相等時間間隔的時間序列影像集，本研究以10 d內(nèi)所有高質(zhì)量影像的NDVI最大值作為該10 d的NDVI值[18,27]。如果影像受云、雪或其他因素的影響，在某些地區(qū)無法獲得10 d內(nèi)高質(zhì)量影像的NDVI值[8,17,22]，本研究則采用相鄰像元線性插值算法來填補缺失值，以構建完整的NDVI時間序列數(shù)據(jù)集[28]。

即使經(jīng)過嚴格的預處理，NDVI時間序列數(shù)據(jù)集中仍殘留由云、大氣和二向性反射引起的噪聲[29]。因此，本研究利用S-G濾波器對NDVI時間序列曲線進行平滑，以消除時間序列中的噪聲[30]。

S-G濾波器于1964年被首次提出[31]，是一種在給定滑動窗口寬度內(nèi)基于局部多項式最小二乘擬合的濾波方法，該濾波器最大的特點為在濾除噪聲的同時可以保持曲線的變化趨勢。使用S-G濾波器對NDVI時間序列進行平滑如式（2）所示。

式中Y為平滑后的NDVI值；Y為滑動窗口內(nèi)平滑前Y前后第個點的NDVI值；C為濾波系數(shù)，即第期Y在滑動窗口中的權重；為滑動窗口的數(shù)據(jù)個數(shù)2+1，其中為滑動窗口寬度。

1.3.2 確定時間窗口

為確定提取越冬作物的最佳時間窗口，本研究選取具有代表性的越冬作物、春播作物、落葉森林、常綠森林、水體和建筑用地的樣本點創(chuàng)建NDVI時間序列曲線。本研究設定滑動窗口為90 d，使用S-G濾波對NDVI時間序列曲線逐像元進行平滑。由于水體和建筑用地的NDVI曲線在整個研究期（2017年9月1日至2018年6月30日）保持穩(wěn)定，水體的NDVI<0，建筑用地的NDVI<0.4，與植被曲線差異較大，因此主要對越冬作物和其他植被的NDVI曲線差異進行研究。

由平滑后越冬作物、春播作物、落葉森林和常綠森林的NDVI曲線可知（圖4），越冬作物在10月中下旬開始播種，此時NDVI值最?。∟DVI<0.3），11月份初出苗后，NDVI值逐步增大，1月份停止生長進入越冬期后，NDVI值保持相對不變，2月份越冬作物開始返青，3月份進入生長旺盛期，NDVI值恢復增長，并于成熟期達到峰值（NDVI>0.6），此后NDVI值逐步下降，到5月底至6月初完成收割，NDVI變?yōu)樽钚≈担∟DVI<0.3）；春播作物在2月中下旬播種前，9月份至第二年2月中下旬NDVI值較低（NDVI<0.3），播種后NDVI值逐漸增大，6月份下旬達到最大值（NDVI>0.6）；落葉森林在10月份之前NDVI值較高（NDVI>0.6），10月份后NDVI值逐漸降低，并于12月份達到NDVI最低值（NDVI<0.3），3月份上旬返青后NDVI逐漸增大（NDVI>0.6）；常綠森林在整個時期均保持較高NDVI值（NDVI>0.6）。

由以上分析得出，越冬作物的NDVI曲線與其他植被的NDVI曲線之間存在明顯差異，即越冬作物播種后生長旺盛，春播作物的耕地撂荒、落葉森林枯萎（NDVI值較低）。越冬作物在播種和收獲時，春播作物、落葉森林和常綠森林則處于生長旺盛期（NDVI值較高）。因此，本研究選擇越冬作物生長旺盛期、播種期和收獲期來區(qū)分越冬作物和其他植被。根據(jù)越冬作物生命周期，確立越冬作物生長旺盛期為2018年03月20日至2018年04月20日，越冬作物播種期為2017年10月11日至2017年11月10日，越冬作物收獲期為2018年05月20日至2018年06月30日。

1.3.3 分類參數(shù)

對于越冬作物生長旺盛期，在Google Earth Engine云平臺使用Max函數(shù)逐像元求NDVI最大值（NDVImax），獲得NDVImax圖像。對于越冬作物播種期和收獲期，在Google Earth Engine云平臺使用Median函數(shù)和Min函數(shù)逐像元分別求解2個時期的NDVI中位值（NDVImedian）和NDVI最小值（NDVImin），獲得NDVImedian和NDVImin圖像。本研究使用NDVImedian是由于在部分云陰影影響區(qū)，常綠森林的NDVI值較低，會與越冬作物混淆，故使用NDVImedian提高越冬作物分類閾值，避免這類型常綠森林錯誤地分成越冬作物，提高了分類精度。至此，獲得了NDVImax、NDVImedian和NDVImin3個分類參數(shù)。

1.3.4 精度驗證方法

本研究基于158個驗證樣本，對本分類模型進行精度評價。首先，以每個驗證樣本點為中心，生成100 m×100 m的矩形地面樣方，共獲取樣方158個（包括越冬作物65個，其他地類93個）；其次，獲取地面樣方空間分辨率為1 m的Google Earth衛(wèi)星影像，對每個地面樣方的影像進行目視解譯以獲得地表參考數(shù)據(jù)；最后，將地表參考數(shù)據(jù)與本研究分類結(jié)果構建混淆矩陣，來評價本研究分類算法的精度。評價指標包括用戶精度（User Accuracy，UA，%）、生產(chǎn)者精度（Producer Accuracy，PA，%）、總精度（Overall Accuracy，OA，%）和Kappa系數(shù)，其計算如式（3）～式（6）所示

式中n為混淆矩陣中第行列的值，n為混淆矩陣中第行的和，為混淆矩陣中第列的和，為驗證樣本總數(shù)，為混淆矩陣行列數(shù)。

2 結(jié)果與分析

2.1 分類模型構建結(jié)果

為了構建分類模型，本研究基于463個訓練樣本（186個越冬作物樣本點、65個春播作物樣本點、57個落葉森林樣本點、50個常綠森林樣本點、90個建筑用地樣本點和15個水體樣本點），繪制了NDVImax、NDVImedian和NDVImin3個分類參數(shù)的箱型圖如圖5所示，通過分析越冬作物與其他地類的特征參數(shù)差異來構建決策樹分類模型。具體操作如下：

1）由于選取的樣本數(shù)無法完全代表研究區(qū)內(nèi)所有越冬作物的物候特征，為了確保能提取出所有越冬作物，本研究將閾值的范圍向外擴大了10%。由圖5a可知，越冬作物的NDVImax的范圍在(0.36,0.82)之間，擴大后的范圍為(0.33,0.91)，而水體和建筑用地的NDVImax上限小于0.33，因此將越冬作物NDVImax的閾值設置為0.33可將水體和建筑用地的像元排除。

2）去除了水體和建筑用地像元的干擾后，下一步將區(qū)分越冬作物和其他植被。由圖5b可知，越冬作物NDVImedian的范圍為(0.19,0.45)，擴大后的范圍為(0.17,0.50)，有75%其他植被的像元值在0.50之上。因此將NDVImedian的閾值設置為0.50去除75%的其他植被像元。由圖5a可知，落葉森林和春播作物的NDVImax范圍分別在[0.24,0.44]和[0.18,0.43]之間，擴大范圍后這2類植被的NDVImax最高不超過0.48，超過0.48的落葉森林和春播作物像元將被去除。由于越冬作物NDVImax的范圍在(0.36,0.82)之間，因此分為NDVImax≤0.48和NDVImax>0.48兩種情況來處理：1）當研究區(qū)內(nèi)像元的NDVImax≤0.48時，可以通過設置NDVImin閾值來去除常綠森林像元。由圖5c可知，常綠森林的NDVImin范圍在0.44以上，此時越冬作物像元的NDVImin范圍在(-0.11,0.16)之間，擴大后為(-0.12,0.17)，因此設置判別條件-0.12< NDVImin<0.17來去除該條件下常綠森林像元；2）當研究區(qū)內(nèi)像元的NDVImax>0.48時，設置判別條件-0.12 < NDVImin<0.5NDVImax去除剩余的常綠森林像元。

3）本研究還選取高程模型中的坡度數(shù)據(jù)來提升精確度，越冬作物設置的坡度閾值為10°。

基于上述分析，構建越冬作物決策樹分類模型如圖6所示。使用決策樹模型第一層“坡度小于10°”判定條件去除非耕地像元，使用決策樹模型第二層“NDVImax>0.33”判定條件去除水體和建筑用地像元，使用決策樹模型第三層“NDVImedian<0.5”判定條件去除落葉森林和春播作物像元。在決策樹模型的最后一層，若像元滿足“NDVImax≤0.48”和“-0.12< NDVImin<0.17”條件，或滿足“NDVImax>0.48”和“-0.12< NDVImin< 0.5NDVImax”條件，則該像元被判定為越冬作物，否則為其他地類。

a. NDVI最大值 a. NDVImaxb. NDVI中位值 b. NDVImedianc. NDVI最小值 c. NDVImin

注：NDVImax為越冬作物生長旺盛期（2018年3月20日至2018年4月20日）NDVI最大值；NDVImedian和NDVImin分別為越冬作物播種期（2017年10月11日至2017年11月10日）和收獲期（2018年5月20日至2018年6月30日）中NDVI中位值和最小值。

Note: NDVImaxis the maximum value of NDVI during the period of peak growth stage (from 2018-03-20 to 2018-04-20); NDVImedianand NDVIminare the median and minimum values of NDVI during the period of sowing stage (from 2017-10-11 to 2017-11-10) and harvest stage (from 2018-05-20 to 2018-06-30),respectively.

圖5 各地類樣本的分類參數(shù)箱型圖

Fig.5 Classification parameter box plots of samples of various land types

2.2 越冬作物分類結(jié)果

淮河流域越冬作物的總面積為8.762×106hm2，其中淮河流域內(nèi)山東省、河南省、安徽省和江蘇省的作物種植面積分別為1.299×106、2.895×106、2.431×106和2.126×106hm2（圖7）。在空間分布方面，32°～35°N是越冬作物的主要分布區(qū)域，占淮河流域越冬作物總面積的89%；113°～120°E是越冬作物的主要分布區(qū)域，占淮河流域越冬作物總面積的97%；淮河流域東北部、西南部和西部越冬作物分布較少。調(diào)查發(fā)現(xiàn)東北部和西部多為山脈（圖1），地勢高，土質(zhì)較差，且農(nóng)田多為梯田，面積較小、收割較困難、不易大規(guī)模種植作物，大多數(shù)農(nóng)田只在夏季種植一季作物（例如玉米、生姜、芋頭、番薯）后便撂荒，冬季不再種植。研究區(qū)西南地區(qū)（例如信陽市）氣溫高降水足，水塘密集，灌溉充分，多種植水稻。水稻在9月末至10月初收割，水稻收割期通常雨水較多，導致稻田較濕，影響了越冬作物的播種，故大部分區(qū)域種植一季，僅有地勢較高區(qū)域種植兩季，淮河流域內(nèi)部空白處主要為城市區(qū)、果園區(qū)（如碭山梨園區(qū)）、湖泊和河流（洪澤湖和微山湖）等。

2.3 精度驗證結(jié)果

基于1.3.4節(jié)獲取的地表參考數(shù)據(jù)，包括528個越冬作物像元（47.521 hm2）和819個其他像元（73.714 hm2），計算混淆矩陣?；煜仃嚨脑蕉魑锓诸惤Y(jié)果精度驗證如表1所示，本研究構建的模型對研究區(qū)越冬作物分類的用戶精度為92.6%，生產(chǎn)者精度為97.0%，總精度95.8%，分類精度較高；Kappa系數(shù)為0.912，說明分類結(jié)果與地表參考數(shù)據(jù)二者間的一致性較強。

將基于目視解譯后的地表參考數(shù)據(jù)與本研究模型的分類結(jié)果進行直觀對比如圖8所示，所選中的Google Earth影像區(qū)域是包含典型地物的矩形區(qū)域（圖 8a），該區(qū)域內(nèi)除越冬作物外，還包括房屋村莊、魚塘、河流和森林等非越冬作物。由圖8b和圖8c對比可知，本研究模型的分類結(jié)果與目視解譯后的地表參考數(shù)據(jù)相似度較高，表明本研究構建的模型以較高的精度提取出了研究區(qū)的越冬作物。分類誤差（包括錯誤分類和遺漏分類）如圖8d所示，越冬作物錯誤分類現(xiàn)象主要集中在耕地路網(wǎng)上，這是由于農(nóng)村道路寬度一般在30 m以內(nèi)，本研究構建的30 m空間分辨率數(shù)據(jù)集對這類線狀地物識別能力較差；越冬作物遺漏分類則主要集中在村落建筑物周圍，村落邊界的不規(guī)則導致鑲嵌在建筑物周邊的小型越冬作物種植區(qū)域無法正確識別。因此，本研究構建的越冬作物提取模型能夠滿足大區(qū)域分類精度的需求。

2019年河南省統(tǒng)計年鑒[34]中縣級單元（即本研究選出的河南省53個縣市）越冬作物種植面積與本研究模型提取的越冬作物種植面積之間的相關性分析如圖9所示。由圖9可知，縣級年鑒統(tǒng)計的種植面積與本研究提取的種植面積具有顯著的線性關系（2=0.86，<0.01），表明本研究模型的分類結(jié)果與統(tǒng)計數(shù)據(jù)有較高的一致性。其中縣級年鑒統(tǒng)計數(shù)據(jù)中驗證區(qū)域的越冬作物總種植面積為3.030×106hm2，本研究提取的越冬作物在驗證區(qū)域的總種植面積為2.804×106hm2，精度為92.6%。

表1 越冬作物分類結(jié)果精度驗證

3 結(jié) 論

本研究基于Google Earth Engine云平臺，融合Landsat-7/8和Sentinel-2A/B影像獲取高時空分辨率影像集，結(jié)合越冬作物獨特的物候期，構建了一種簡單、有效和高精度的越冬作物識別模型，在像元尺度上繪制淮河流域越冬作物種植面積圖，主要的結(jié)論如下：

1）融合多個傳感器影像能夠獲取高時空分辨率數(shù)據(jù)集，克服當前大面積作物制圖所用數(shù)據(jù)源時空分辨率的缺陷（例如MODIS空間分辨率低、Landsat時間分辨率低等）。本研究結(jié)合了所有可使用的高質(zhì)量Landsat-7/8和Sentinel-2 A/B衛(wèi)星影像（合計4 982景影像），極大地增加了每個像元上可用的高質(zhì)量影像觀測數(shù)，更好地捕獲了越冬作物的物候。

2）通過10 d合成、插值和平滑算法，重構了時間序列歸一化植被指數(shù)（Normalized Difference Vegetation Index，NDVI）數(shù)據(jù)集，揭示了越冬作物和其他植被的NDVI曲線差異，確定了越冬作物關鍵物候期：生長旺盛期（2018年3月20日至4月20日）和播種期（2017年10月11日至11月10日）和收獲期（2018年5月20日至6月30日），成功構建了融合越冬作物獨特物候特征的分類模型。

3）2018年淮河流域越冬作物種植面積為8.762×106hm2，其精度滿足大區(qū)域作物制圖的精度要求。

淮河流域處于中國東部季風區(qū)，受云的影響，很難保證每個像元10 d有一幅高質(zhì)量衛(wèi)星影像，本研究采用相鄰像元線性插值算法來填補缺失值，但是在拐點處缺失影像時，插值并不能完全反映作物生長的真實位置，存在一定的不確定性。Sentinel-1A/B影像不受天氣的影像，能夠估算作物的物候期，未來可結(jié)合Sentinel-1 A/B衛(wèi)星的C波段合成孔徑雷達（Synthetic Aperture Radar，SAR）數(shù)據(jù)進行制圖。本研究構建的30 m空間分辨率數(shù)據(jù)集難以提取研究區(qū)地塊面積小于0.09 hm2的農(nóng)田，未來可基于高分1/2/3/5/6等衛(wèi)星獲取10 m或更高空間分辨率的作物分布圖。此外，云平臺的高性能計算和并行處理能力可用于處理大尺度的地理空間數(shù)據(jù)集，是未來遙感大數(shù)據(jù)處理研究的趨勢。

致謝：本研究得到了河南大別山森林生態(tài)系統(tǒng)國家野外科學觀測研究站和國家地球系統(tǒng)科學數(shù)據(jù)中心-黃河中下游分中心（http://henu.geodata.cn）的數(shù)據(jù)支持，在此表示謝意！

[1] Franch B, Vermote E F, Becker-Reshef I, et al. Improving the timeliness of winter wheat production forecast in the United States of America, Ukraine and China using MODIS data and NCAR growing degree day information[J]. Remote Sensing of Environment, 2015, 161: 131-148.

[2] Shao Y, Campbell J B, Taff G N, et al. An analysis of cropland mask choice and ancillary data for annual corn yield forecasting using MODIS data[J]. International Journal of Applied Earth Observation and Geoinformation, 2015, 38: 78-87.

[3] 王利民，劉佳，楊福剛，等. 基于 GF-1 衛(wèi)星遙感的冬小麥面積早期識別[J]. 農(nóng)業(yè)工程學報，2015，31(11)：194-201.

Wang Limin, Liu Jia, Yang Fugang, et al. Early recognition of winter wheat area based on GF-1 satellite[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(11): 194-201. (in Chinese with English abstract)

[4] Hunt M L, Blackburn G A, Carrasco L, et al. High resolutions wheat yield mapping using Sentinel-2[J/OL]. Remote Sensing of Environment, 2019, 233, [2019-09-07], https: //doi. org/10. 1016/j. rse. 2019. 111410.

[5] Jiao X F, Kovacs J, Shang J L, et al. Object-oriented crop mapping and monitoring using multi-temporal polarimetric RADARSAT-2 data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 96: 38-46.

[6] Jin Z N, Azzari G, Lobell D. Improving the accuracy of satellite-based high-resolution yield estimation: A test of multiple scalable approaches[J]. Agricultural and Forest Meteorology, 2017, 247: 207-220.

[7] 解毅，張永清，荀蘭，等. 基于多源遙感數(shù)據(jù)融合和LSTM算法的作物分類研究[J]. 農(nóng)業(yè)工程學報，2019，35(15)：129-137.

Jie Yi, Zhang Yongqing, Xun Lan, et al. Crop classification based on multi-source remote sensing data fusion and LSTM algorithm[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(15): 129-137. (in Chinese with English abstract)

[8] Liu L, Xiao X M, Qin Y W, et al. Mapping cropping intensity in China using time series Landsat and Sentinel-2 images and Google Earth Engine[J/OL]. Remote Sensing of Environment, 2020, 239, [2019-12-20], https: //doi. org/10. 1016/j. rse. 2019. 111624.

[9] Liu C, Zhang Q, Tao S Q, et al. A new framework to map fine resolution cropping intensity across the globe: Algorithm, validation, and implication[J/OL]. Remote Sensing of Environment. 2020, 251, [2020-09-08], https: //doi. org/10. 1016/j. rse. 2020. 112095.

[10] Tian H F, Huang N, Niu Z, et al. Mapping winter crops in China with multi-source satellite imagery and phenology-based algorithm[J]. Remote Sensing, 2019, 11(7): 820-843.

[11] 周亮，任建強，吳尚蓉，等. 基于卷積神經(jīng)網(wǎng)絡的中國北方冬小麥遙感估產(chǎn)[J]. 農(nóng)業(yè)工程學報，2019，35(15)：119-128.

Zhou Liang, Ren Jianqiang, Wu Shangrong, et al. Remote sensing estimation on yield of winter wheat in North China based on convolutional neural network[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(15): 119-128. (in Chinese with English abstract)

[12] 張錦水，趙光政，洪友堂，等. 基于像元物候曲線匹配的生長季內(nèi)河北省冬小麥空間分布識別[J]. 農(nóng)業(yè)工程學報， 2020，36(23)：193-200.

Zhang Jinshui, Zhao Guangzheng, Hong Youtang, et al. Spatial extraction of winter wheat in Hebei in growing season using pixel-wise phenological curve[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(23): 193-200. (in Chinese with English abstract)

[13] Ren S J, Guo B, Wu X, et al. Winter wheat planted area monitoring and yield modeling using MODIS data in the Huang-Huai-Hai Plain, China[J/OL]. Computers and Electronics in Agriculture, 2021, 182, [2021-02-07], https:// doi. org/10. 1016/j. compag. 2021. 106049.

[14] 黃健熙，趙劍橋，汪雪淼，等. 基于遙感和積溫的冬小麥生育期提取方法[J]. 農(nóng)業(yè)機械學報， 2019，50(2)：176-183.

Huang Jianxi, Zhao Jianqiao, Wang Xuemiao, et al. Extraction method of growth stages of winter wheat based on accumulated temperature and remote sensing data[J]. Transactions of the Chinese Society for Agricultural Machinery, 2019, 50(2): 176-183. (in Chinese with English abstract)

[15] Phalke A, ?zdo?an M, Thenkabail P, et al. Mapping croplands of Europe, Middle East, Russia, and Central Asia using Landsat, Random Forest, and Google Earth Engine[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2020, 167: 104-122.

[16] Wang S, Tommaso S D, Deines J M, et al. Mapping twenty years of corn and soybean across the US Midwest using the Landsat archive[J/OL]. Scientific Data, 2020, 7, [2020-08-13], https: //doi. org/10. 1038/s41597-020-00646-4.

[17] Pan L, Xia H M, Zhao X Y, et al. Mapping winter crops using a phenology algorithm, time-series Sentinel-2 and Landsat-7/8 images, and Google Earth Engine[J/OL]. Remote Sensing, 2021, 13, [2021-06-21], https: //doi. org/10. 3390/rs13132510.

[18] Griffiths P, Nendel C, Hostert P. Intra-annual reflectance composites from Sentinel-2 and Landsat for national-scale crop and land cover mapping[J]. Remote Sensing of Environment, 2019, 220: 135-151.

[19] Li R Y, Xu M Q, Chen Z Y, et al. Phenology-based classification of crop species and rotation types using fused MODIS and Landsat data: The comparison of a random-forest-based model and a decision-rule-based model[J/OL]. Soil and Tillage Research. 2021, 206, [2020-10-07], https: //doi. org/10. 1016/j. still. 2020. 104838.

[20] Lenco D, Interdonato R, Gaetano R, et al. Combining Sentinel-1 and Sentinel-2 Satellite Image Time Series for land cover mapping via a multi-source deep learning architecture[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2019, 158: 11-22.

[21] 王情，劉雪華，岳天祥. 淮河流域糧食生產(chǎn)潛力空間格局研究[J]. 生態(tài)經(jīng)濟，2014，30(7)：24-27.

Wang Qing, Liu Xuehua, Yue Tianxiang. Study on the spatial patterns of food production potential in the Huaihe River Basin[J]. Ecological Economy, 2014, 30(7): 24-27. (in Chinese with English abstract)

[22] Pan L, Xia H M, Yang J, et al. Mapping cropping intensity in Huaihe Basin using phenology algorithm, all Sentinel-2 and Landsat images in Google Earth Engine[J/OL]. International Journal of Applied Earth Observation and Geoinformation, 2021, 102, [2021-05-20], https: //doi. org/10. 1016/j. jag. 2021. 102376.

[23] Wulder M, Masek J, Cohen W, et al. Opening the archive: How free data has enabled the science and monitoring promise of Landsat[J]. Remote Sensing of Environment, 2012, 122: 2-10.

[24] Torres R, Snoeij P, Geudtner D, et al. GMES Sentinel-1 mission[J]. Remote Sensing of Environment, 2012, 120: 9-24.

[25] Foga S, Scaramuzza P, Guo S, et al. Cloud detection algorithm comparison and validation for operational Landsat data products[J]. Remote Sensing of Environment, 2017, 194: 379-390.

[26] 王世炎，夏雨春，等. 河南統(tǒng)計年鑒2019[M]. 北京：中國統(tǒng)計出版社，2019.

[27] Holben B. Characteristics of maximum-value composite images from temporal AVHRR data[J]. International Journal of Remote Sensing, 1986, 7(11): 1417-1434.

[28] Chen Y L, Song X D, Wang S S, et al. Impacts of spatial heterogeneity on crop area mapping in Canada using MODIS data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2016, 119: 451-461.

[29] Cihlar J, Ly H, Li Z Q, et al. Multitemporal, multichannel AVHRR data sets for land biosphere studies-Artifacts and corrections[J]. Remote Sensing of Environment, 1997, 60(1): 35-57.

[30] Wang J, Xiao X M, Liu L, et al. Mapping sugarcane plantation dynamics in Guangxi, China, by time series Sentinel-1, Sentinel-2 and Landsat images[J/OL]. Remote Sensing of Environment. 2020, 247, [2020-06-09], https: //doi. org/10. 1016/j. rse. 2020. 111951.

[31] Savitzky A, Golay M. Smoothing and differentiation of data by simplified least squares procedures[J]. Analytical Chemistry, 1964, 36(8): 1627-1639.

Mapping of the winter crop planting areas in Huaihe River Basin based on Google Earth Engine

Pan Li1, Xia Haoming1,2※, Wang Ruimeng1, Niu Wenhui1, Tian Haifeng1, Qin Yaochen1,2

(1.,-,475001,; 2.(),,475001,)

The winter crop has been one of the important crop types in China. Accurate and timely spatio-temporal distribution of planting area directly determines the grain output and economy, as well as the national food security. Taking the Huaihe River Basin as an example, this study aims to extract the planting areas of winter crops according to the phenology period using the Google Earth Engine cloud platform and the fusion of Landsat-7/8 and Sentinel-2A/B images. Firstly, a dataset of time-series images was constructed with a spatial resolution of 30 m. A CFMask algorithm was selected to preprocess the images, thereby calculating the Normalized Difference Vegetation Index (NDVI). More importantly, the maximum NDVI of all high-quality images within 10 days was used to obtain time-series data with equal time intervals. The linear interpolation was utilized to fill the pixels without high-quality images. Savitzky-Golay (S-G) filtering (a second-order filter with a moving window of 9 observations) was adopted to smooth the NDVI time series for the removal of noise. As such, a smoothed NDVI time series was obtained with a 10-day interval. Secondly, the peak growth, sowing, and harvest periods were determined to select sample points of winter crops with different spatial distributions, according to the NDVI time series. Subsequently, the winter crops were sowed in mid-late October, when the NDVI values were the lowest. The NDVI values gradually increased, after the emergence of seedlings in early November. The crops stopped growing in January during the overwintering period, where the NDVI stayed the same over the whole period. Furthermore, the NDVI resumed growing and gradually reached the peak growth period, when the winter crops turned green in February. After that, the NDVI reached the peak at the heading stage, and then gradually decreased. Correspondingly, the NDVI dropped to the bottom, when the harvest was over from the end of May to June. According to these characteristics in the process of winter crops growth, the peak growth period was determined from March 20, 2018, to April 20, 2018, the sowing period was determined from October 11, 2017, to November 10, 2018, and the harvest period was determined from May 20, 2018, to June 30, 2018. Particularly, the maximum NDVI was achieved in the peak growth period and the minimum and median of NDVI in the sowing and harvest period. Finally, the classification model of a decision tree was constructed, according to the NDVI boxplots of winter crops and non-winter crops at different time periods. The planting area map of winter crops was also generated for the Huaihe River Basin. The results showed that the planting area of winter crops was 8.762×106hm2in the Huaihe River Basin in 2018. Specifically, the user accuracy was 0.926, the producer accuracy was 0.970, the total accuracy was 0.958, and the Kappa coefficient was 0.912. Consequently, the large-scale planting area of winter crops was extracted accurately for the decision-making in similar areas.

remote sensing; crops; mapping; planting area; Google Earth Engine; Huaihe River Basin

潘力，夏浩銘，王瑞萌，等. 基于Google Earth Engine的淮河流域越冬作物種植面積制圖[J]. 農(nóng)業(yè)工程學報，2021，37(18)：211-218.doi：10.11975/j.issn.1002-6819.2021.18.025 http://www.tcsae.org

Pan Li, Xia Haoming, Wang Ruimeng, et al. Mapping of the winter crop planting areas in Huaihe River Basin based on Google Earth Engine[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(18): 211-218. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2021.18.025 http://www.tcsae.org

2020-11-10

2021-08-30

河南省科技攻關計劃項目（212102310019）；國家自然科學基金項目（41701433）；黃河文明省部共建協(xié)同創(chuàng)新中心重大項目（2020M19）；河南省自然科學基金資助項目（202300410531）

潘力，研究方向為農(nóng)業(yè)遙感。Email: panli970611@henu.edu.cn

夏浩銘，博士，副教授，研究方向為定量遙感及其綜合應用。Email: xiahm@vip.henu.edu.cn

10.11975/j.issn.1002-6819.2021.18.025

TP79

A

1002-6819(2021)-18-0211-08