基于時序Sentinel-2的河北地下水漏斗區作物單雙季種植遙感監測

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

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

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

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

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

1 試驗部分

1.1 研究區概況

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

基于中國氣象數據網(http://data.cma.cn)和河北省統計局提供的中國農作物生長發育數據集,歸納總結得出研究區范圍內1年1熟作物和1年2熟作物主要物候信息見圖2。

1.2 數據獲取與預處理

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

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

1.3 試驗方法

1.3.1 技術路線 本研究技術路線見圖3。

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

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

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

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

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

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

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

隨著冬季農作物的生長發育,其色彩飽和度(saturation,S)逐漸增加,而其他地物并不會有明顯的S變化或者隨著氣溫的降低,S會略微降低,故可以采用S去除土壤等易混淆地物。根據實測數據,采集農作物與非農作物12月的S和11月的S并作差,得到的即為飽和度S的變化值,農作物S差值主要集中在0.05以上,而非農作物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系數和F1分數。F1分數同時兼顧了分類模型的精確率和召回率。其公式如下:

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

2 結果與分析

2.1 作物周年光譜變化曲線

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

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

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

2.2 不同時序數據組合遙感監測精度

不同時序遙感數據組合分類場景作物單雙季種植監測結果見圖6,其F1分數指數見表1。

表1 不同時序數據組合場景F1分數精度評價

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

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

2.3 不同分類方法遙感監測精度

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

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

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

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

2.4 作物單雙季種植監測結果

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

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

2.5 地下水漏斗區作物種植調整建議

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

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

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

3 結論與討論

3.1 結論

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

3.2 討論

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