基于Landsat-8和Sentinel-2時間序列合成影像的山區甘蔗種植區提取

康倩文　徐偉恒　王雷光　洪澤湖　劉運

關鍵詞：甘蔗；物候；時間序列影像；影像合成；谷歌地球引擎

中圖分類號：S566.1 文獻標識碼：A

甘蔗是全球第一大糖料作物和第二大生物能源作物，中國是世界上最主要的甘蔗種植國家之一[1]。從甘蔗中提取的蔗糖是中國食品行業中不可替代的重要原料，而云南省作為我國第二大蔗糖生產基地[2]，蔗糖產業是云南省的高原特色產業，更是扶貧產業[3]。云南的甘蔗主要種植于丘陵山區旱坡地，準確識別甘蔗種植區分布范圍，對有關部門制定相應的土地政策，為甘蔗產量和效益評估提供有效參考，對當地農業可持續發展、保障食糖供給安全都具有重要現實意義。

農作物的遙感識別方法主要包括兩大類：單時相遙感影像識別和多時相遙感變化監測[4]。多時相遙感影像相較單時相影像更能反映植被的季節變化以及物候特征，有效增加農作物的識別精度，目前已經被廣泛應用于農作物的識別與監測中[5]。如邊增淦等[6]利用2018 年18 景GF-1 WFV影像構建出的歸一化植被指數時間序列并采用分層策略逐步提取出黑河流域中游地區的作物種植結構。劉吉凱等[7]選用多時相GF_1 WFV 高分辨率遙感數據，采用多時相迭代法構建甘蔗提取決策樹模型提取出廣西江州區的甘蔗種植區，總體精度（overall accuracy， OA）達到95.26%。大量研究表明，聯合多源遙感影像可以增加影像數量，增加對地觀測有效像元次數[8]，有效彌補單種數據源數據量少的局限性[9]，是實現多云霧地區、大范圍農作物種植信息提取的主要手段[10]。如譚深等[11]融合Landsat-8 與Sentinel-2 光學數據、Sentinel-1 微波SAR 數據和其他輔助數據，解決了南方地區多云天氣導致影像質量不佳的問題，準確、高效地提取了海南省的水稻種植范圍，OA為93.2% 。ZHENG 等[12] 基于Landsat-7/8 、Sentinel-1/2 合成的時間序列影像利用TWDTW法評估作物的物候相似性，并識別出中國的甘蔗種植區，SR 數據的OA 為93.47%。

不同作物具有不同的物候特征，植被物候能夠反映植被生長規律，國內外學者通過提取植被物候特征，進行甘蔗等農作物的識別與分類。ALEXANDRE 等[13]利用多時相MODIS-EVI 數據通過甘蔗與其他作物的增強型植被指數物候特征的差異實現巴西圣保羅州農作物的分類。張東東等[14]通過選擇特定時相的HJ-1 CCD 影像，利用甘蔗生長特點，采用決策樹分類模型區分出中國南方地區的甘蔗與其他地物，OA 達到93.2%，Kappa 系數為0.81。WANG 等[15]基于Sentinel-1、Sentinel-2 和Landsat 影像構建的時間序列數據，采用基于像元和物候學的方法繪制了2018 年廣西省甘蔗分布圖， 甘蔗的制圖精度（produceraccuracy， PA）和用戶精度（user accuracy， UA）分別為88%和96%。

綜上所述，目前學者們已經利用多時相、多數據源并結合物候特征的方法對甘蔗種植區進行了提取，但是缺乏專門針對山區甘蔗提取的研究。由于山區地物類型復雜，耕地碎片化現象嚴重，耕地地塊細小狹窄且結構復雜[16]，導致耕地信息難以快速、精準獲取。為明確山區甘蔗種植區分布的情況和規律，本研究綜合利用多源、多時相遙感數據與植被物候的特點，以云南省玉溪市新平縣為研究區，采用谷歌地球引擎（Google EarthEngine， GEE）云計算平臺，基于Landsat-8 和Sentinel-2 長時間序列合成影像數據，依據甘蔗與不同地物之間的物候特征、地形特征，借助物候參數，探究適用于山區甘蔗種植區提取的方法。旨在開發一種適合提取山區甘蔗種植區的有效算法，為當地農業部門制定合理的土地利用政策提供參考。

1 材料與方法

1.1 研究區概況

新平縣位于云南省玉溪市西部，地形復雜，多山地、丘陵，耕地空間破碎化，種植結構復雜，是云南省甘蔗主要生產縣之一，也是云南元江干熱河谷的典型山區農業縣[17] 。新平縣地處23°38′15″～24°26′05″N，101°16′30″～102°16′50″E，全縣面積約為4223 km²（圖1）。研究區地形以山地為主，地處哀牢山中段東麓，縣境山區面積達4139.6 km²，占縣域面積的98%；地勢西北高、東南低，海拔為418～3143 m。新平縣年均降雨800 mm 左右，降雨量少，蒸發量大，氣候干熱，全年無霜，有利于蔗糖的積累[18]，適宜栽種甘蔗。蔗糖產業是新平縣的傳統產業，是當地農民增收、企業增效的主要來源之一。

1.2 甘蔗提取技術路線

首先，將Landsat-8 OLI 和Sentinel-2 MSI 影像數據經過去云、光譜波段校準后計算遙感指數得到遙感指數時間序列，并將影像進行合成、插值、平滑，構建出以10 d 為間隔的高質量無云、無空洞的時間序列合成影像；其次，利用時間序列合成影像提取耕地圖層，并基于耕地圖層，利用甘蔗的物候特征、地形特征提取甘蔗，并進行精度驗證（圖2）。

1.3 數據源及預處理

1.3.1 遙感影像及預處理 本研究基于GEE 平臺調用了2019 年10 月1 日到2021 年7 月1 日的覆蓋研究區所有的Landsat-8 OLI 與Sentinel-2MSI 影像共計1022 張。所用Landsat-8 OLI 影像、Sentinel-2 MSI 影像分別來自美國地質調查局提供的Landsat-8 一級表面反射率數據集和Sentinel-2 的Level-2A 級地球表觀反射率數據集（表1），2 種數據集均已經過幾何校正、輻射定標、大氣校正， 同時利用Landsat-8 OLI 的QA_PIXEL 波段和Sentinel-2 MSI 的QA 波段分別對數據集進行云和云陰影的掩膜。研究區內像元的總觀測次數的空間分布見圖3A，研究區內像元的良好觀測次數即掩膜云和云陰影后的像元的有效觀測次數的空間分布見圖3B，統計結果顯示，所有像元在2019 年10 月1 日到2021 年7 月1 日期間至少有35 次良好的觀測，圖3A 與圖3B 中右下角觀測次數較高的原因是該區域恰好位于Landsat-8 OLI 行列號為129/43、129/44、130/43的影像的重疊區域以及Sentinel-2 MSI 行列號為N0208_R104 、N0209_R104 、N0300_R104 、N0208_R061、N0209_R061、N0300_R061 的影像的重疊區域。

本研究采用ZHANG 等[19]提供的校正系數對Landsat-8 和Sentinel-2 影像進行光譜波段校準，以消除不同衛星傳感器之間的差異。計算每張影像的光譜指數： 歸一化植被指數（ normalizeddifference vegetation index， NDVI）[20]、陸表水分指數（land surface water index， LSWI）[21]、增強型植被指數（enhanced vegetation index， EVI）[22]和修正后歸一化差異水體指數（modified normalizeddifference water index， MNDWI）[23]，計算公式如式（1）～（4）所示。

1.3.2 數字高程模型 數字高程模型（digitalelevation model， DEM）可生成坡度、坡向、剖面圖等，輔助地形地貌分析。本研究所采用的DEM開源數據來源于美國國家航空航天局提供的第三級產品， 空間分辨率為30 m ， 產品ID 為“USGS/SRTMGL1_003”。由于新平縣多山地的地形地貌，因此本研究考慮引入地形因子來輔助分析甘蔗的地形特征。

1.3.3 野外調查數據 為確定新平縣甘蔗的物候特征，本團隊于2022 年4 月進行實地調研。樣本集由實地調查定位出在2020 年作物的種植類型未發生改變的樣本和谷歌地球結合目視解譯的解譯樣本組成，包括常綠植被、水體、不透水層、甘蔗與其他農作物。甘蔗與其他農作物的樣本組成如表2 所示。

1.4 時間序列影像合成

影像合成壓縮了數據的時間維度，減少了數據冗余，合成后的時間序列影像十分規則有序[24]，有利于研究地物的物候特征。最大值合成法可避免云雨天氣對遙感影像計算的NDVI 值存在低值噪聲的影響[25]，因此，本研究以10 d 為影像合成的時間窗口，以最大值合成法生成NDVI 時間序列，以均值合成法生成LSWI、EVI、MNDWI 時間序列，并將合成后的影像全部重采樣為10 m。采用線性插值法[25]填補影像由于去云和去云影造成的空洞，利用移動窗口大小為9，濾波器階數為2 的Savitzky-Golay 濾波器來平滑NDVI 時間序列，以消除數據噪聲影響[26]。基于此，本研究構建了新平縣2019 年10 月1 日至2021 年7 月1 日的時間序列合成影像數據集。

1.5 甘蔗種植區提取算法

實地調研顯示，新平縣的主要土地類型大致分為4 類：常綠植被、水體、不透水層和耕地。借鑒前人的理論基礎通過遙感指數時間序列對4類地物進行提取，識別算法如表3 所示。其中，常綠植被在時間序列中被定義為全年存在超過90%的觀測值的LSWI 大于0 且EVI 大于0.2；水體在時間序列中被定義為全年存在至少5%的觀測值的MNDWI 大于NDVI 且NDVI 小于0.1；不透水層在時間序列中被定義為全年超過90%的觀測值的LSWI 小于0.2。耕地圖層通過掩膜常綠植被、水體、不透水層的圖層范圍得到，并通過在耕地圖層中研究甘蔗與其他作物的物候差異來實現對甘蔗種植區的提取。

植物物候現象是環境條件季節和年際變化最直觀、最敏感的生物指示器[30]。植被指數時間序列能夠較好地反映植被的代謝強度及其季節性變化及年際變化。通過研究分析甘蔗與研究區內主要農作物（水稻、玉米、香蕉）的NDVI 時間序列曲線（圖4），發現新平縣甘蔗大多在3 月末至4 月初播種，此時甘蔗的NDVI 值最低，9 月甘蔗的NDVI 值上升到0.7 以上，11—12 月甘蔗成熟，12 月至次年4 月為甘蔗收割期，甘蔗的NDVI 值開始下降；研究區內水稻大多為雙季稻，在1 年內有2 個生長周期，與甘蔗明顯不同；玉米比甘蔗的種植時間晚而收獲時間早；香蕉與甘蔗的NDVI 曲線類似，但是甘蔗收割后表現出與裸土相似的NDVI 值，與香蕉明顯不同。甘蔗的生長周期較長，一般在8 個月至1 年以上[31]，并且甘蔗是一種高生物量作物，因此甘蔗有較長一段時間的NDVI 值處于較高水平。

為了定量研究甘蔗與其他農作物的物候特征差異，我們引入物候參數來說明。植被指數時間序列被廣泛應用于植被物候參數提取[32]。利用物候參數可以較好地進行農作物分類研究[33]。本研究根據NDVI 時間序列曲線提取的物候參數如圖5 所示，物候參數的說明如表4 所示，其中基值（base of season， BOS）為生長季階段NDVI 曲線左邊最小值（a 點）與右邊最小值（f 點）的均值，振幅（amplitude， AMP）為生長季階段NDVI 曲線最大值與基值的差值，上升時間（rise time， Tr）表示生長季階段NDVI 從振幅的10%（b 點）上升到振幅的90%（c 點）所用的時間，下降時間（fall time， Tf）表示生長季階段NDVI 從振幅的90%（d 點）下降到振幅的10%（e 點）所用的時間，上積分（above integral of season， AIOS）表示生長季階段NDVI 曲線與基值之間的積分，下積分（below integral of season， BIOS）表示生長季階段基值與坐標0 軸之間的積分[34]。本研究選擇物候參數上升時間、下降時間、上積分、下積分作為甘蔗提取的物候特征，并繪制甘蔗與其他農作物的訓練樣本在物候參數上升時間、下降時間、上積分、下積分上的直方圖分布（圖6）。本研究基于90%的置信區間獲取甘蔗訓練樣本的物候參數范圍，并依據此范圍進行研究區內甘蔗種植區的提取。物候參數閾值分布范圍為：上升時間大于30 d 且小于190 d，下降時間大于70 d 且小于200 d，上積分大于90 且小于164，下積分大于26 且小于82。

由于新平縣地形以山地為主，因此本研究考慮引入海拔、坡度因子作為山區甘蔗提取的輔助特征。據研究報道，云南省甘蔗海拔分布范圍為422～1300 m[35-37]，坡度達到25°以上。實地調研顯示， 新平縣甘蔗種植區大致分布在海拔480～1330 m，坡度30°以下。因此綜合考慮研究報道以及實地調研結果調整研究區內提取甘蔗種植區的海拔及坡度因子范圍為海拔為400～1400 m，坡度在30°以下。

綜上，本研究提出的甘蔗種植區提取算法為：30

2 結果與分析

2.1 主要土地類型提取結果

本研究提取的新平縣常綠植被、水體、不透水層和耕地如圖7A 所示。圖7B 為圖7A 中標記為1 區域的細節放大圖，綠色區域為新平縣常綠植被提取結果的局部示意圖；圖7C 為圖7A 中標記為2 區域的細節放大圖，藍色區域為新平縣水體提取結果的局部示意圖；圖7D 為圖7A 中標記為3 區域的細節放大圖，粉色區域為新平縣不透水層提取結果的局部示意圖；圖7E 為圖7A 中標記為4 區域的細節放大圖，黃色區域為新平縣耕地提取結果的局部示意圖。由圖7A 可知，新平縣耕地主要集中在新平縣中部及北部。

2.2 甘蔗提取結果及精度評價

為了驗證本研究算法對研究區甘蔗的提取精度，依據實地調研情況在Google Earth 中隨機選取274 個甘蔗樣本與470 個其他農作物樣本作為驗證樣本，基于混淆矩陣構建OA、Kappa 系數、PA 和UA 進行山區甘蔗提取精度評估。新平縣甘蔗提取的OA 達97.07%，Kappa 系數為0.83，從分類精度評價指標來看，提取OA 較高，整體分類精度滿足空間分析與實際應用需求。甘蔗提取的PA 與UA 分別達到88.86%和80.57%，其他農作物PA 與UA 分別達到97.88%和98.89%，提取面積為7705 hm²，提取結果較好。

本研究提取的2020 年新平縣的甘蔗種植區分布圖見圖8A，其中圖8B 和圖8C 分別為圖8A中標記為1、2 區域的細節放大圖（透明度60%的淺綠色區域為識別的甘蔗）。由圖8A 可知，2020年新平縣甘蔗種植區大多分布于新平縣中南部及北部，集中于漠沙鎮、戛灑鎮、老廠鄉，其他地區甘蔗種植區零星分布。從甘蔗種植區分布地形地貌來看，2020 年新平縣甘蔗種植區多分布于新平縣北部及中部偏南，分布于地勢較平坦、海拔相對較低的地帶，由于新平縣地勢西北高、東南低，甘蔗種植區也呈現東南多于西北的現象。進一步分析研究區內甘蔗種植區的空間分布情況可知，新平縣甘蔗地塊分散且不規則，地塊邊界模糊。該分布圖可以作為開展甘蔗種植區長勢、災害等遙感監測的基礎數據。從鄉鎮尺度來看，新平縣甘蔗分布面積較多的區域是漠沙鎮、戛灑鎮、老廠鄉，其中漠沙鎮甘蔗分布最多（圖9）。漠沙鎮、戛灑鎮、老廠鄉的甘蔗種植面積分別為2786、1458、1332 hm²，分別占新平縣甘蔗種植面積的36.16%、18.92%、17.29%。新平縣甘蔗分布面積較少的區域分別是古城街道、桂山街道、建興鄉，古城街道甘蔗分布最少，古城街道、桂山街道、建興鄉的甘蔗種植面積分別為0.87、54.00 、13.70 hm²，分別占新平縣甘蔗種植面積的0.01%、0.70%、0.18%。分析發現新平縣各鄉鎮甘蔗提取分布趨勢情況與實地調研結果一致。

2.3 物候參數敏感性分析

對本研究中用到的上升時間、下降時間、上積分、下積分4 個物候參數進行敏感性分析（圖10）。可以發現，選取上升時間、下降時間、上積分、下積分這4 個物候參數作為提取山區甘蔗種植區的物候特征時，甘蔗提取算法中甘蔗的UA為88.85%達到最高，PA 為80.57%，此時甘蔗的錯分誤差最小，漏分誤差最大。當甘蔗提取算法中物候參數分別設為下降時間、上積分、下積分（Tf+AIOS+BIOS），上升時間、上積分、下積分（Tr+AIOS+BIOS），上升時間、下降時間、上積分（Tr+Tf+AIOS），上升時間、下降時間、下積分（Tr+Tf+BIOS）這4 種情況時，甘蔗提取算法中甘蔗的UA 分別下降為88.50%、79.37%、75.35%、67.20%，PA 分別上升為82.42%、84.40%、83.52%、85.00%。物候參數敏感性分析表明，選取上升時間、下降時間、上積分、下積分4 個物候參數參與甘蔗提取算法時甘蔗提取效果最好，4個物候參數在提取甘蔗時都是重要的物候參數，都是甘蔗制圖算法的重要組成部分。本研究甘蔗制圖算法中所用物候參數對于耕地中甘蔗種植區的提取有較好的提取能力，有助于甘蔗與其他農作物的物候特征分析。

3 討論

本研究提出了適用于山區復雜地形，基于Landsat-8 OLI 與Sentinel-2 MSI 合成的時間序列遙感影像數據的甘蔗種植區的提取方法。以典型山區甘蔗種植區新平縣為研究區，生成一幅10 m空間分辨率的2020 年新平縣甘蔗種植區分布圖。本研究中甘蔗種植區提取的OA 為97.07%，Kappa系數為0.83， 提取精度較為良好， 表明基于Landsat-8 OLI 與Sentinel-2 MSI 合成的時間序列遙感影像數據結合光譜指數特征、物候特征與地形特征可以實現山區甘蔗種植區的提取并繪制出甘蔗分布圖，精度滿足基本需求。提取的新平縣甘蔗種植區面積為7705 hm²，其中漠沙鎮、戛灑鎮、老廠鄉的甘蔗種植最多，總共占新平縣甘蔗種植區的72.37%，古城街道、桂山街道和建興鄉的甘蔗種植最少，總共占新平縣甘蔗種植區的0.89%，可見鄉鎮間甘蔗種植區分布差異顯著。本研究提出的甘蔗種植區提取算法的UA和PA 不是很高，其原因是新平縣耕地種植結構復雜，地塊破碎，中等空間分辨率的影像像元內可能包含多種地類；本研究選擇可免費獲取的Landsat-8 和Sentinel-2 影像在多云天氣的新平縣很難保證每個像元在10 d 內有一幅高質量的衛星影像，經過插值后的像元值在反映作物生長的真實情況時存在一定的不確定性，因此未來可利用不受天氣影響的Sentinel-1 A/B 影像來輔助進行甘蔗種植區的提取。

新平縣是典型的山區甘蔗種植區，在利用物候特征提取山區甘蔗時，多源、多時相遙感數據的融合可極大提高遙感影像的空間與時間分辨率，盡可能避免山區多云雨、復雜的地形地貌因素導致山區遙感影像質量不高[38]的問題，為山區甘蔗提取與精細化制圖提供一種有效參考。本研究基于NDVI 時間序列數據，借助甘蔗與常綠植被、水體、不透水層、其他農作物在光譜指數特征、物候特征、地形特征上的差異，采用上升時間、下降時間、上積分、下積分4 個物候參數以及海拔、坡度因子確定提取甘蔗的最佳閾值，相較張東東等[14]學者選擇用多時相影像，單一光譜特征，利用純凈甘蔗訓練樣本進行決策樹分類的方法來進行甘蔗種植區的提取，本研究不僅增加了研究區內像元良好觀測次數，克服了山區等多云雨地區影像質量不高的問題，并且可以更充分地利用甘蔗的生長特征及物候特征，對甘蔗樣本長勢情況的依賴相對較小。本研究利用4 個物候參數確定了甘蔗提取的最佳物候參數閾值，上升時間、下降時間、上積分、下積分4 個物候參數在提取甘蔗時都是重要的物候參數，為未來其他地區的山區甘蔗種植區提取提供重要的物候特征參考。

研究區內地貌復雜，耕地破碎度較高，由于研究區海拔及氣候等因素的影響，耕地在休耕期會出現雜草叢生，作物收割殘留枝葉等現象，在影像上導致光譜混淆，容易出現“異物同譜”現象。因此，使得在甘蔗種植區提取過程中難免出現錯提或漏提現象。今后考慮將使用更高分辨率的遙感數據或運用更小的時間分辨率組成的時間序列數據來進一步提高提取的精度，建立普適性更高的山區甘蔗種植區提取算法。本研究的提取精度雖然滿足縣域的山區甘蔗提取精度，適用于作為當地農業部門的決策建議，但是甘蔗提取的用戶精度和制圖精度還是有待提高。希望之后考慮通過引入甘蔗的紋理特征等，提高山區甘蔗的提取精度，實現山區甘蔗的精細化提取。