基于時序Sentinel-2的河北地下水漏斗區(qū)作物單雙季種植遙感監(jiān)測

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

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

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

隨著遙感大數(shù)據(jù)時代的來臨,高空間分辨率的遙感數(shù)據(jù)全球陸續(xù)免費開放獲取,為農(nóng)業(yè)領域帶來了新的選擇[5-6]。遙感提取作物種植范圍研究發(fā)展有近50年時間,目前農(nóng)作物種植制度遙感提取的方法主要有多時序遙感數(shù)據(jù)曲線分析和遙感分類2類。多時序遙感數(shù)據(jù)曲線分析法主要依靠長時間序列光譜分析結合物候信息進行不同熟制農(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é)點和分類特征構建分類場景,篩選農(nóng)作物關鍵時間節(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)地物種類眾多,分布零散,氣候類型為溫帶大陸性季風氣候,適合多種作物生長。研究區(qū)主要農(nóng)作物為小麥、玉米、棉花,有少量大豆,種植制度僅包含1年1熟和1年2熟[13]。

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

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

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

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

1.3 試驗方法

1.3.1 技術路線 本研究技術路線見圖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ù)據(jù),結合周年衛(wèi)星影像數(shù)據(jù),對比1年1熟作物和1年2熟作物的光譜曲線,挑選出6個關鍵時間節(jié)點,進行作物單雙季種植遙感監(jiān)測最佳時相組合試驗。關鍵時間節(jié)點分別為:3個峰值,K1:12月,K2:4—5月,K3:7—8月;1個谷值,K4:6月;2個作物不同趨勢時間節(jié)點,K5:11月,K6:3月。

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

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

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

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

根據(jù)實測數(shù)據(jù),采集冬季作物的NDVI和H,構建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熟作物的閾值設置為NDVI>0.33,H處于100～220之間,12月和11月S差值大于0。

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

1.4 精度評價

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

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

2 結果與分析

2.1 作物周年光譜變化曲線

本研究基于GEE遙感云平臺,以月為單位組合拼接30 d內(nèi)高質(zhì)量影像,結合野外調(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)測的關鍵時間節(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)測結果見圖6,其F1分數(shù)指數(shù)見表1。

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

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

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

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

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

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

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

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

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

基于NDVI、EVI、NDWI、GCVI、STI、NDTI等6種植被指數(shù)和哨兵2號數(shù)據(jù)的B3、B4、B8、B11這4個波段共10個特征,結合場景6中6個關鍵時間節(jié)點共8個月的時序組合數(shù)據(jù)進行分類監(jiān)測,獲得河北平原地下水漏斗區(qū)1年1季作物和1年2季作物分布圖,對研究區(qū)內(nèi)2塊田塊進行分類細節(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)象,地塊較為零碎,分類效果最差;最大似然法分類結果尚可,區(qū)分出了區(qū)域內(nèi)絕大多數(shù)地塊,但是地塊零碎,連續(xù)性較差,誤分情況較多。

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

為更加直觀地比較不同區(qū)域單雙季作物種植分布情況,結合河北省縣市圖,制作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田塊尺度單雙季作物分類結果、河北平原漏斗區(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 結論與討論

3.1 結論

本研究基于2020年10月至2021年9月的周年哨兵2號遙感數(shù)據(jù),通過隨機森林法監(jiān)測河北平原地下水漏斗區(qū)作物單雙季種植情況,得到以下結論:(1)1年2熟作物和1年1熟作物周年光譜曲線之間的差異主要有6個關鍵時間節(jié)點,且使用6個關鍵時間節(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ù)且考慮云層覆蓋問題,對關鍵時間節(jié)點的影像抓取有一定的難度,也會對區(qū)分精度產(chǎn)生一定的影響。(3)由于研究區(qū)面積廣、東西南北有較大跨度,不同的氣候條件會導致各個地區(qū)作物生長速率、耕作時間存在一定的差異,影響提取精度,后續(xù)可以考慮針對同緯度大范圍農(nóng)場進行調(diào)查分析以提升精度。