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

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

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

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

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

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

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

0 引 言

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

1.2 數(shù)據(jù)來(lái)源

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

Landsat-7/8衛(wèi)星攜帶增強(qiáng)型專題制圖儀（Enhanced Thematic Mapper，ETM+）和陸地成像儀（Operational Land Image ，OLI），空間分辨率為30 m，重訪周期為15 d[23]。Landsat-7 ETM+和Landsat-8 OLI影像分別對(duì)應(yīng)“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對(duì)應(yīng)“COPERNICUS / S2”。本研究以越冬作物的一個(gè)完整物候期為例，選取時(shí)間范圍為2017年9月1日至2018年6月30日，合計(jì)4 982景影像。

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

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

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

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

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

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

1.3 研究方法

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

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

基于NDVI時(shí)間序列曲線識(shí)別作物是常見(jiàn)且有效的方法，NDVI的計(jì)算如（1）所示。

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

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

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

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

1.3.2 確定時(shí)間窗口

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

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

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

1.3.3 分類參數(shù)

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

1.3.4 精度驗(yàn)證方法

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

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

2 結(jié)果與分析

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

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

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

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

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

基于上述分析，構(gòu)建越冬作物決策樹(shù)分類模型如圖6所示。使用決策樹(shù)模型第一層“坡度小于10°”判定條件去除非耕地像元，使用決策樹(shù)模型第二層“NDVImax>0.33”判定條件去除水體和建筑用地像元，使用決策樹(shù)模型第三層“NDVImedian<0.5”判定條件去除落葉森林和春播作物像元。在決策樹(shù)模型的最后一層，若像元滿足“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為越冬作物生長(zhǎng)旺盛期（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），地勢(shì)高，土質(zhì)較差，且農(nóng)田多為梯田，面積較小、收割較困難、不易大規(guī)模種植作物，大多數(shù)農(nóng)田只在夏季種植一季作物（例如玉米、生姜、芋頭、番薯）后便撂荒，冬季不再種植。研究區(qū)西南地區(qū)（例如信陽(yáng)市）氣溫高降水足，水塘密集，灌溉充分，多種植水稻。水稻在9月末至10月初收割，水稻收割期通常雨水較多，導(dǎo)致稻田較濕，影響了越冬作物的播種，故大部分區(qū)域種植一季，僅有地勢(shì)較高區(qū)域種植兩季，淮河流域內(nèi)部空白處主要為城市區(qū)、果園區(qū)（如碭山梨園區(qū)）、湖泊和河流（洪澤湖和微山湖）等。

2.3 精度驗(yàn)證結(jié)果

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

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

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

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

3 結(jié) 論

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

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

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

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

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

致謝：本研究得到了河南大別山森林生態(tài)系統(tǒng)國(guó)家野外科學(xué)觀測(cè)研究站和國(guó)家地球系統(tǒng)科學(xué)數(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)星遙感的冬小麥面積早期識(shí)別[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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è)工程學(xué)報(bào)，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] 周亮，任建強(qiáng)，吳尚蓉，等. 基于卷積神經(jīng)網(wǎng)絡(luò)的中國(guó)北方冬小麥遙感估產(chǎn)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，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] 張錦水，趙光政，洪友堂，等. 基于像元物候曲線匹配的生長(zhǎng)季內(nèi)河北省冬小麥空間分布識(shí)別[J]. 農(nóng)業(yè)工程學(xué)報(bào)， 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è)機(jī)械學(xué)報(bào)， 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)濟(jì)，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)計(jì)年鑒2019[M]. 北京：中國(guó)統(tǒng)計(jì)出版社，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è)工程學(xué)報(bào)，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

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

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

夏浩銘，博士，副教授，研究方向?yàn)槎窟b感及其綜合應(yīng)用。Email: xiahm@vip.henu.edu.cn

10.11975/j.issn.1002-6819.2021.18.025

TP79

A

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