鄭翔宇, 李存軍, 劉 玉, 孟浩然, 楊鐵利, 盧 闖, 張方娟

(1.北京市農(nóng)林科學(xué)院信息技術(shù)研究中心,北京 100097; 2.遼寧科技大學(xué)土木工程學(xué)院,遼寧鞍山 114051;3.農(nóng)業(yè)部農(nóng)村農(nóng)業(yè)遙感機理與定量遙感重點實驗室,北京 100097)

作為重要糧食主產(chǎn)區(qū)的河北平原種植制度主要以1年2熟種植方式為主,1年1熟種植方式為輔。1年2熟作物對水資源的剛性需求,決定了小麥生育期70%～80%的水需要通過灌溉來滿足,導(dǎo)致河北平原近年來每年地下水以1.0～1.5 m的速度下降[1-2]。研究表明,春玉米生長期與華北平原降水的時空分布耦合度最好,可基本實現(xiàn)生育期最小灌溉或零灌溉[3]。準確獲取大區(qū)域田塊尺度單雙季作物種植信息對于種植結(jié)構(gòu)調(diào)整、農(nóng)業(yè)水資源保護和合理利用有重要價值[4]。

隨著遙感大數(shù)據(jù)時代的來臨,高空間分辨率的遙感數(shù)據(jù)全球陸續(xù)免費開放獲取,為農(nóng)業(yè)領(lǐng)域帶來了新的選擇[5-6]。遙感提取作物種植范圍研究發(fā)展有近50年時間,目前農(nóng)作物種植制度遙感提取的方法主要有多時序遙感數(shù)據(jù)曲線分析和遙感分類2類。多時序遙感數(shù)據(jù)曲線分析法主要依靠長時間序列光譜分析結(jié)合物候信息進行不同熟制農(nóng)作物空間分布情況獲取。黃翀利用中分辨率成像光譜儀(MODIS)數(shù)據(jù)獲取水稻歸一化植被指數(shù)(NDVI)年時間序列曲線,提取柬埔寨水稻熟制信息[7];Xiao等基于MODIS數(shù)據(jù)多源信息融合方法獲得了南亞和東南亞單雙季水稻農(nóng)業(yè)制圖,該方法多采用中、低分辨率的多光譜影像,在提取大尺度區(qū)域性農(nóng)作物方面具有明顯優(yōu)勢[8]。遙感分類方法可以通過閾值或遙感分類器獲取不同種植制度農(nóng)作物空間分布情況。王雙等基于SPOT-5衛(wèi)星數(shù)據(jù),采用最大似然法分類的方式,對河北省三河市單季玉米和棉花進行提取[9];Son等基于哨兵2號(Sentinel-2)衛(wèi)星數(shù)據(jù)采用遙感分類方法對我國臺灣地區(qū)單雙季水稻種植面積進行提取[10],該方法在提取精度方面具有明顯優(yōu)勢,但由于大氣、云層等因素的影響,較難獲得長時序的農(nóng)作物空間分布情況且計算量較大。然而,河北平原地區(qū)作物遙感檢測過去多采用250 m MODIS長時間序列分析方法[11],空間分辨率低且混合像元嚴重,缺乏大范圍田塊尺度作物種植制度監(jiān)測信息;如何選擇合理的遙感時間序列組合減少遙感數(shù)據(jù)量以提高分類效率有待進一步研究[12]。

本研究利用10 m空間分辨率的Sentinel-2影像,提取并分析河北平原地下水漏斗區(qū)作物全年NDVI曲線,篩選分類時間節(jié)點和分類特征構(gòu)建分類場景,篩選農(nóng)作物關(guān)鍵時間節(jié)點組成最優(yōu)時序特征組合,制作大范圍田塊尺度農(nóng)作物單雙季種植分類圖,并分析和提出地下水漏斗區(qū)作物種植調(diào)整建議。

1 試驗部分

1.1 研究區(qū)概況

研究區(qū)為河北平原地下水漏斗區(qū)。利用國家30 m數(shù)字高程模型,計算河北省南部地區(qū)的坡度值,基于坡度和高程將河北省南部平原地區(qū)提取出來作為研究區(qū)。研究區(qū)位于河北省南部(圖1),地處36°40′～39°10′N,114°30′～117°00′E,研究區(qū)面積為66 020 km2。研究區(qū)內(nèi)地物種類眾多,分布零散,氣候類型為溫帶大陸性季風(fēng)氣候,適合多種作物生長。研究區(qū)主要農(nóng)作物為小麥、玉米、棉花,有少量大豆,種植制度僅包含1年1熟和1年2熟[13]。

基于中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn)和河北省統(tǒng)計局提供的中國農(nóng)作物生長發(fā)育數(shù)據(jù)集,歸納總結(jié)得出研究區(qū)范圍內(nèi)1年1熟作物和1年2熟作物主要物候信息見圖2。

1.2 數(shù)據(jù)獲取與預(yù)處理

本研究基于谷歌地球引擎(google earth engine,GEE)平臺獲取Sentinel-2 Level-2A級產(chǎn)品。選取2020年10月至2021年9月的影像數(shù)據(jù),進行作物單雙季種植遙感監(jiān)測。Sentinel-2影像經(jīng)過拼接、裁剪、去云、波段計算和融合等預(yù)處理工作,獲得Sentinel-2的耕地時間序列數(shù)據(jù)集。

2020年11月對研究區(qū)域進行實地調(diào)查和走訪,共獲取實地調(diào)查樣本點650個,其中包括1年1熟作物地塊樣點309個、1年2熟作物地塊樣點341個,按照7 ∶3的比例將其分為訓(xùn)練樣本和驗證樣本,然后將待分類場景與訓(xùn)練樣本輸入到分類器中[14]。

1.3 試驗方法

1.3.1 技術(shù)路線 本研究技術(shù)路線見圖3。

1.3.2 光譜和植被指數(shù) 根據(jù)研究區(qū)作物種植情況,本研究選用NDVI、增強植被指數(shù) (enhanced vegetation index,EVI)、歸一化水指數(shù) (normalized difference water index,NDWI)、葉綠素植被指數(shù) (green chlorophyll vegetation index,GCVI)、土壤耕作指數(shù) (soil tillage index,STI)、差異耕作指數(shù)(normalized difference tillage index,NDTI)6種植被指數(shù)和哨兵2號數(shù)據(jù)的B3、B4、B8、B11這4個波段數(shù)據(jù)作為特征集,用于作物單雙季種植遙感監(jiān)測。

1.3.3 分類場景設(shè)計 本研究通過地面實測數(shù)據(jù),結(jié)合周年衛(wèi)星影像數(shù)據(jù),對比1年1熟作物和1年2熟作物的光譜曲線,挑選出6個關(guān)鍵時間節(jié)點,進行作物單雙季種植遙感監(jiān)測最佳時相組合試驗。關(guān)鍵時間節(jié)點分別為:3個峰值,K1:12月,K2:4—5月,K3:7—8月;1個谷值,K4:6月;2個作物不同趨勢時間節(jié)點,K5:11月,K6:3月。

本研究設(shè)計9種不同時間序列特征組合的分類場景(圖3)。

1.3.4 分類方法 本研究采用隨機森林算法進行作物單雙季種植遙感監(jiān)測,并與傳統(tǒng)的閾值法和最大似然法對比。

隨機森林法是基于決策樹的自學(xué)習(xí)集成式機器學(xué)習(xí)算法,運算速度快、模型穩(wěn)健、泛化能力強,是農(nóng)作物遙感制圖的常用技術(shù)手段。

閾值法采用多時相遙感數(shù)據(jù)和HSV(色調(diào)、飽和度、明度)變化提取1年2熟作物范圍[15-16],采用NDVI閾值法提取夏季作物范圍并減去1年2熟范圍得到1年1熟作物分布范圍。將哨兵2號數(shù)據(jù)的B4、B8、B11這3個波段的數(shù)據(jù)融合,并轉(zhuǎn)化為HSV(色調(diào)、飽和度、明度)數(shù)據(jù)以分離出色調(diào)(hue,H)和明度(value,V)的數(shù)據(jù),通過小麥和非小麥地物在H-NDVI空間上的差異,采用HSV閾值法初步識別冬小麥疑似區(qū)域。

根據(jù)實測數(shù)據(jù),采集冬季作物的NDVI和H,構(gòu)建H-NDVI圖(圖4),可以看出研究區(qū)內(nèi)冬季作物主要分布在NDVI大于0.33且H處于 100～220之間的區(qū)域,采用此閾值對其進行初步提取。

隨著冬季農(nóng)作物的生長發(fā)育,其色彩飽和度(saturation,S)逐漸增加,而其他地物并不會有明顯的S變化或者隨著氣溫的降低,S會略微降低,故可以采用S去除土壤等易混淆地物。根據(jù)實測數(shù)據(jù),采集農(nóng)作物與非農(nóng)作物12月的S和11月的S并作差,得到的即為飽和度S的變化值,農(nóng)作物S差值主要集中在0.05以上,而非農(nóng)作物S差值主要集中在0以下。綜上所述,提取1年2熟作物的閾值設(shè)置為NDVI>0.33,H處于100～220之間,12月和11月S差值大于0。

1年1熟作物種植范圍閾值提取,通過NDVI對7—8月耕地掩膜后的影像進行提取,閾值設(shè)置為NDVI>0.45,提取結(jié)果減去1年2熟作物范圍得到1年1季作物范圍。

1.4 精度評價

精度評價采用構(gòu)建混淆矩陣的方法,計算分類結(jié)果的制圖精度(PA)、用戶精度(UA)、總體精度(OA),Kappa系數(shù)和F1分數(shù)。F1分數(shù)同時兼顧了分類模型的精確率和召回率。其公式如下:

其中:N表示訓(xùn)練樣本的個數(shù);n表示分類類別數(shù);i表示某一具體類別;Xjj代表混淆矩陣某一類地物正確分類的個數(shù);X+j代表該類別真實參考總數(shù)(混淆矩陣該地物X對應(yīng)列的總和);Xii代表正確分到該類的像元總數(shù)(混淆矩陣對角線值);Xi+代表整個影像分類為該地物的像元個數(shù)(混淆矩陣該地物X對應(yīng)行的總和);X+i代表整個影像正確分類為該地物的像元個數(shù);Xij代表所有像元個數(shù);OA指在驗證樣本集上分類器預(yù)測正確的概率;F1分數(shù)由用戶精度和制圖精度計算得到,是用戶精度和制圖精度調(diào)和均值,用于表征模型輸出的好壞,取值范圍位于 0～1之間,其值越接近1表示模型輸出越好;Kappa系數(shù)是檢測驗證樣本與預(yù)測結(jié)果吻合度的指標,其值的大小通常位于0～1之間,越接近1表示吻合度越高[17]。

2 結(jié)果與分析

2.1 作物周年光譜變化曲線

本研究基于GEE遙感云平臺,以月為單位組合拼接30 d內(nèi)高質(zhì)量影像,結(jié)合野外調(diào)查,制作的1年1熟作物和1年2熟作物NDVI變化曲線見圖5。

1年1熟作物(如春玉米、棉花)在4—5月開始陸續(xù)種植,NDVI開始升高,在7—8月NDVI達到頂峰后于9—10月陸續(xù)收割,NDVI迅速降低,在10月中下旬收割完成,周年NDVI光譜曲線呈單峰趨勢。1年2熟作物在秋季末期開始種植,并在10—11月份迅速生長,在冬季初期(12月)形成小峰值,由于冬季雨、雪、氣溫降低等氣候因素的影響,NDVI在1—2月有所降低,在初春時期開始繼續(xù)生長發(fā)育直到4月份達到第2個小高峰。1年2熟作物在5—6月成熟和收割,并在6—7月重新種植夏季作物,故在6月呈現(xiàn)裸地情況,顯現(xiàn)出NDVI光譜曲線的深谷,又在7—8月迅速生長發(fā)育,呈現(xiàn)出NDVI最高峰值現(xiàn)象。故1年2熟作物NDVI周年光譜曲線呈現(xiàn)1年3峰情況,與1年1熟作物有較大差別,時序遙感信息可以幫助有效區(qū)分2種種植制度。

對比2種光譜曲線,遙感監(jiān)測的關(guān)鍵時間節(jié)點分為峰值、谷值和不同趨勢時間節(jié)點3種情況,第1個峰值在12月,第2個峰值在4—5月,第3個峰值在7—8月,谷值在6月;趨勢時間節(jié)點是指1年1熟作物和1年2熟作物NDVI增減趨勢不同的時間節(jié)點,主要有3—4月時1年1熟作物還沒開始種植,NDVI無明顯變化,而1年2熟作物NDVI快速增加的時間節(jié)點;11—12月1年1熟作物已經(jīng)收割完成,NDVI無明顯變化,而1年2熟作物開始生長發(fā)育,NDVI開始增加。

2.2 不同時序數(shù)據(jù)組合遙感監(jiān)測精度

不同時序遙感數(shù)據(jù)組合分類場景作物單雙季種植監(jiān)測結(jié)果見圖6,其F1分數(shù)指數(shù)見表1。

表1 不同時序數(shù)據(jù)組合場景F1分數(shù)精度評價

由圖6、表1可知,區(qū)分1年1熟作物和1年2熟作物的最佳分類場景為場景6,即所有關(guān)鍵時間節(jié)點(8個月)組成的特征時序,該場景對1年1熟作物和1年2熟作物分類的F1分數(shù)指數(shù)都是最高的。

分析場景1～場景7(表1)可以發(fā)現(xiàn),在關(guān)鍵時間節(jié)點峰值場景中,峰值時間序列對區(qū)分1年1熟作物和1年2熟作物具有較好的效果,3個峰值關(guān)鍵時間節(jié)點的精度要高于2個峰值關(guān)鍵時間節(jié)點(場景1～場景2)。對比谷值和2個趨勢不同時間節(jié)點,其中11—12月的趨勢不同關(guān)鍵時間節(jié)點影響較大,對區(qū)分精度的影響略高于3—4月的關(guān)鍵時間節(jié)點和6月的谷值時間節(jié)點(場景3～場景5)。對比所有關(guān)鍵時間節(jié)點時序融合影像和周年時序融合影像的區(qū)分精度(場景6～場景7),僅使用關(guān)鍵時間節(jié)點的時序遙感數(shù)據(jù)精度略微高于周年遙感數(shù)據(jù)精度,說明關(guān)鍵時間節(jié)點的時序遙感數(shù)據(jù)能明顯地將2種種植制度進行區(qū)分,可能是由于研究區(qū)范圍非常大,同一作物存在物候差異,且周年遙感影像數(shù)據(jù)無法做到完全無云情況,故會對區(qū)分精度產(chǎn)生一定的影響,導(dǎo)致精度下降。

2.3 不同分類方法遙感監(jiān)測精度

利用傳統(tǒng)的閾值提取方法、最大似然法與隨機森林方法的作物單雙季種植遙感監(jiān)測結(jié)果見圖7,F1得分精度見表2。其中最大似然法采用所有關(guān)鍵時間節(jié)點融合影像,訓(xùn)練集和驗證集同隨機森林分類場景一致。

從圖7可以看出,無論是針對1年1熟作物提取還是1年2熟作物提取,隨機森林算法的精度優(yōu)于閾值法和最大似然法,最大似然法的區(qū)分精度次之,閾值法的區(qū)分精度最低。

通過對比不同分類方法的F1分數(shù)(表2),隨機森林算法的分類精度最優(yōu),閾值法和最大似然法相對較低,這種情況是由于不同方法的分類原理所導(dǎo)致的。閾值法是通過對影像光譜分析和數(shù)據(jù)分析制定合理的閾值,區(qū)分不同類別,閾值的設(shè)置和光譜數(shù)據(jù)的選擇都會對精度產(chǎn)生較大的影響;由于研究區(qū)跨度較大,研究區(qū)范圍內(nèi)地物種類繁多且分布不均勻,增加了通過似然度區(qū)分地物的難度,導(dǎo)致最大似然分類法精度較低;而隨機森林算法可以在一定程度上避免一次分類錯誤導(dǎo)致的誤差,且隨機森林分類的多棵決策樹更容易把握光譜值上的細微差異,可以獲取較高精度,對區(qū)分不同耕作方式具有較好的效果。

表2 3種不同分類方法F1分數(shù)精度對比

2.4 作物單雙季種植監(jiān)測結(jié)果

基于NDVI、EVI、NDWI、GCVI、STI、NDTI等6種植被指數(shù)和哨兵2號數(shù)據(jù)的B3、B4、B8、B11這4個波段共10個特征,結(jié)合場景6中6個關(guān)鍵時間節(jié)點共8個月的時序組合數(shù)據(jù)進行分類監(jiān)測,獲得河北平原地下水漏斗區(qū)1年1季作物和1年2季作物分布圖,對研究區(qū)內(nèi)2塊田塊進行分類細節(jié)展示,結(jié)果見圖8。

由圖8可知,1年1熟作物主要分布在研究區(qū)的東北部和南部地區(qū),其中滄州市、廊坊市、邢臺市分布較多,在西部和中部也有少量分布;1年2熟作物主要集中在中部和西部大部分地區(qū),研究區(qū)北部和東部也有少量分布;部分地區(qū)1年1熟作物種植和1年2熟作物種植并存。造成這種情況的原因可能是由于河北平原西部水資源較東部多,1年2熟作物冬季耕作以灌溉為主,需要大量地下水資源,而1年1熟作物分布與河北省地下水漏斗區(qū)域較吻合。由圖8可知,隨機森林算法效果最優(yōu),表現(xiàn)為分類后地塊較為完整,地塊連續(xù)性較好,漏分誤分錯分現(xiàn)象較少,能夠區(qū)分出絕大部分農(nóng)田;該地塊范圍內(nèi)閾值法產(chǎn)生了大量的錯分現(xiàn)象,地塊較為零碎,分類效果最差;最大似然法分類結(jié)果尚可,區(qū)分出了區(qū)域內(nèi)絕大多數(shù)地塊,但是地塊零碎,連續(xù)性較差,誤分情況較多。

2.5 地下水漏斗區(qū)作物種植調(diào)整建議

為更加直觀地比較不同區(qū)域單雙季作物種植分布情況,結(jié)合河北省縣市圖,制作1年2熟作物優(yōu)勢度分布圖(圖9)。優(yōu)勢度是該熟制作物面積與耕地總面積的比值,可以有效描述作物在某一地區(qū)種植面積所占比例[18]。對比地下水缺水區(qū)域大致分布圖(圖10),提出地下水漏斗區(qū)作物種植制度調(diào)整建議,以達到減少地下水消耗的目的。

河北平原地下水漏斗區(qū)1年1熟作物種植面積占所有耕地的45.04%,1年2熟作物種植面積占所有耕地的54.96%。其中,滄州市單季作物占全市耕地的65.18%;石家莊單季作物占全市耕地的31.63%;衡水市單季作物占全市耕地面積的44.21%;保定市單季作物占全市耕地面積的47.60%;廊坊市單季作物占全市耕地面積的88.07%;邯鄲市單季作物占全市耕地面積的41.75%;邢臺市單季作物占全市耕地面積的42.72%。

對比10 m田塊尺度單雙季作物分類結(jié)果、河北平原漏斗區(qū)1年2熟作物優(yōu)勢度圖和河北深層地下水等水埋深分布,可以看出主要地下水漏斗區(qū)作物類型以1年1熟作物為主,1年2熟作物優(yōu)勢度較低,但也有一些漏斗區(qū)作物類型依舊以1年2熟作物為主。圖9中網(wǎng)格線區(qū)域(衡水部分縣市)為地下水漏斗區(qū)的核心地區(qū),但是其1年2熟作物占比較大,需要較大的農(nóng)業(yè)用水,不利于地下水位的保持,建議作為重點管控區(qū)減少1年2熟作物的種植;圖中橫線區(qū)域(衡水、滄州、邢臺部分縣市)為地下水漏斗區(qū)邊緣地區(qū),且其1年2熟作物種植比例略小于網(wǎng)格線區(qū)域,建議作為次要管控區(qū)減少1年2熟作物的種植來保持地下水水位。

3 結(jié)論與討論

3.1 結(jié)論

本研究基于2020年10月至2021年9月的周年哨兵2號遙感數(shù)據(jù),通過隨機森林法監(jiān)測河北平原地下水漏斗區(qū)作物單雙季種植情況,得到以下結(jié)論:(1)1年2熟作物和1年1熟作物周年光譜曲線之間的差異主要有6個關(guān)鍵時間節(jié)點,且使用6個關(guān)鍵時間節(jié)點(共8個月)時序數(shù)據(jù)分類效果最佳,相比全年時序數(shù)據(jù)(12個月)減少了遙感數(shù)據(jù)量、運算時間,同時取得了更高的精度。(2)用3種不同的分類方法開展作物單雙季種植提取,隨機森林方法精度最高,最大似然法精度次之,閾值法精度最低。采用隨機森林算法和最佳時序組合,制作了河北平原地下水漏斗區(qū)10 m空間分辨率田塊尺度的作物單雙季種植分布圖,從圖中可知河北平原地下水漏斗區(qū)1年1熟作物主要分布在東北部和中南部地區(qū),1年2熟作物主要分布在西部和中部地區(qū)。(3)對比地下水漏斗區(qū)1年2熟作物優(yōu)勢度圖與河北平原主要地下水漏斗區(qū)域發(fā)現(xiàn),地下水漏斗區(qū)多為1年1熟作物種植,但衡水大部分地區(qū)、滄州南部和邢臺部分縣市1年2熟作物的種植比例較多,建議減少1年2熟作物的種植面積以控制地下水資源損耗。

3.2 討論

本研究主要針對時間序列特征對區(qū)分1年1熟作物和1年2熟作物進行分類和精度評價,雖然取得了較高的精度,但也存在可優(yōu)化和待解決的問題,有待進一步提升提取精度:(1)冬季農(nóng)作物植被較為稀疏,植被光譜信息與裸地易產(chǎn)生混淆,光譜值較為相似,對提取精度產(chǎn)生了影響,考慮后續(xù)可以針對提取的植被指數(shù)進行調(diào)查分析,篩選更加適合的植被指數(shù)和光譜指數(shù)進行提取。(2)1年2熟作物1年內(nèi)收割和種植時間非常接近,土地空閑時間非常短,采用重訪周期為5 d的哨兵2號數(shù)據(jù)且考慮云層覆蓋問題,對關(guān)鍵時間節(jié)點的影像抓取有一定的難度,也會對區(qū)分精度產(chǎn)生一定的影響。(3)由于研究區(qū)面積廣、東西南北有較大跨度,不同的氣候條件會導(dǎo)致各個地區(qū)作物生長速率、耕作時間存在一定的差異,影響提取精度,后續(xù)可以考慮針對同緯度大范圍農(nóng)場進行調(diào)查分析以提升精度。