時序影像在冬小麥種植區提取中的應用分析

　　摘 要：衛星數據在農業土地信息提取中可以發揮重要作用，而由于農業土地具有顯著的物候差異性，時間序列數據的應用能夠有效提高作物的提取精度。不同作物因為物候期的差異，對時間序列的需求有明顯的差異性。為了確定時間序列數據對冬小麥提取精度的影響，利用 NDVI、EVI 作為構造時序影像的特征，將覆蓋整個冬小麥生育期的所有 Sentinel-2A/B 數據，以三個不同的時間間隔(10 d、15 d 和 30 d)的均值合成，然后使用隨機森林分類器對像素或地塊進行分類，以對比不同時間間隔合成的時序影像在面向不同元素的冬小麥提取時的應用效果和提取效率。結果表明：① 當合成 Sentinel-2A/B 影像的時間間隔縮短時，兩種分類方法的冬小麥的提取精度都會更高，但面向地塊的分類方法能達到更好的提取效果;② 10 d 合成的面向地塊的分類方法獲得了最高的 OA(96.28%)和 Kappa 系數(91.71%);○3 時間序列的合成方法對提取效率影響不大，影響它的主要原因是面向不同元素的分類方法。

　　關鍵詞：圖像分割;時間序列;Google Earth Engine;冬小麥提取;隨機森林

　　石嫻; 明艷芳; 劉春秀; 瞿渝; 隋淞蔓， 無線電工程 發表時間：2021-11-25

　　0 引言

　　冬小麥為頭年秋季播種次年初夏收獲的一種農作物[1]，在全世界范圍內都有著廣闊的種植面積[2]。它的種植面積對經濟和社會發展有重要的影響[3,4]。精確、高效的冬小麥空間分布監測具有十分重要的意義[5]。遙感技術能夠快速實現大面積作物種植區域的連續監測，可以在冬小麥監測中發揮重要作用[6]。

　　農作物的遙感識別方法主要包括兩大類: 單時相遙感影像識別、多時相遙感變化檢測。然而由于不同植被之間具有較為相似的光譜特征，單一的光譜信息難以達到較高的識別精度，而物候期的差異通常被作為農作物識別的重要信息。此外，區域灌溉不均勻、土壤鹽漬化、病蟲害、植被密度差異甚至植物的含水量不同等因素都會導致相同作物在生長過程中產生不同程度的差異[7]，因此很難利用單一時相的遙感影像進行區分。多時相遙感影像能夠充分反映不同作物的物候特性及變化規律，增大光譜特征相似的作物之間的可分離性，可以很好地解決“異物同譜”和“同物異譜”問題，大大提高識別精度[8]。目前隨著高時空分辨率遙感的快速發展，基于時間序列的農作物遙感識別已廣泛應用于作物種植面積提取、農業遙感監測等方面。

　　時間序列遙感作物識別方法主要是利用作物整個生育期固定時間間隔的時序數據獲取農作物特定的時序曲線，然后分析其物候或數學特征并進行農作物識別。楊小喚等繪制并分析了作物的 NDVI 時序曲線，制定了不同作物的提取依據，獲取了北京市的冬小麥、玉米、大豆的種植空間分布圖[9];陳健等以整個作物生育期內 35 個時相的合成地表反射率數據生成 MODIS 數據的 EVI 時序數據，結合作物物候歷和種植結構，構建分類模型，最終達到了 95.7%的整體分類精度[10]。近年來隨著遙感衛星研發技術的不斷進步，云平臺計算能力和計算機圖像處理技術的不斷提高，作物識別逐漸向大空間尺度、長時間序列方向發展，分類算法也在不斷改進。周珂等根據每旬 NDVI 最大值合成 Landsat 8 影像，提取了河南省的冬小麥，其平均分類精度達到了 95.3%，且有效降低了與統計數據的相對誤差[11];You 等人將 3 個波段、7 個光譜指數作為特征，使用 10 天中值合成 Sentinel-2 時間序列影像，利用隨機森林分類器進行分類，首次在 10m 分辨率下制作了中國年度區域作物圖[12]。

　　然而目前大多數的識別方法都是使用某一種時間序列數據，很少有專門的研究通過試驗比較分析不同間隔的時序影像用來做冬小麥提取的適宜性和準確性。根據 Blaes 等人的研究，時間序列影像中過長的時間間隔或關鍵期一景影像的缺失都會對分類精度造成影響[13]。Wardlow 等人表示，每種作物有其獨特的時序光譜曲線[14]，不同的構造方法會造成時序數據的變化，使得作物的提取效果產生差異。由此，時序數據的構造方法對于時序影像的遙感識別分類和基于物候特征的遙感監測產生著至關重要的影響，此類研究也十分有必要開展。另外，由于中等分辨率衛星自身的局限性以及小規模農業種植區域的分類易受到背景因素的干擾，像素尺度下的冬小麥提取效果可能難以滿足要求。面向地塊的分類方法在過往的研究中被證明可以很好地克服上述問題，尤其是在小規模農業盛行的地區要優于面向像素的分類方法[15,16]。但是目前還未有在不同時間序列數據下關于地塊和像素提取方法的對比，也未有關于提取效率方面的評價。

　　綜上所述，冬小麥種植區提取目前主要存在時間序列單一的問題，另外沒有考慮到面向地塊、像素方法與不同時序數據結合的適用性。因此本文利用 GEE 平臺和 Sentinel-2A/B 數據評估了不同時間間隔(10 d、 15 d 和 30 d)合成的時序數據對冬小麥提取精度的影響，明確了更適合提取冬小麥種植區的 Sentinel-2A/B 時序數據。同時對比了使用不同間隔的時序數據時面向地塊和面向像素分類方法的優劣性，為擴大 Sentinel-2A/B 數據在作物分類中的應用提供了依據。

　　1 研究區和數據預處理

　　1.1 研究區概況

　　本文研究區膠州市(36°00′-36°30′N，119°37′-120°12′E)位于山東省青島市，占地面積為 1324km2，總耕地面積 82 萬畝，地貌類型有丘陵、平原、洼地和濱海低地四類。膠州市年無霜期 210d，平均氣溫 12℃，有效積溫 4599.2℃，日照時數 2170.5h，降水量 724.8mm，為典型的暖溫帶半濕潤季風區大陸性氣候，四季分明，雨熱同期，非常適宜小麥、馬鈴薯、花生等農作物的種植。

　　1.2 作物物候期

　　通過實地調查了解到研究區內農作物類型多種多樣，主要種植的農作物有冬小麥、馬鈴薯、花生、玉米等。該地玉米年耕種兩次，冬小麥、馬鈴薯和花生則耕種一次，典型地表類型包括冬小麥-夏玉米雙熟制耕地、馬鈴薯-夏玉米雙熟制耕地、蔬菜-春玉米雙熟制耕地、花生單熟制耕地。冬小麥一般在 10 月初播種， 12 月下旬進入越冬期，翌年 3 月開始返青生長，4 月進入生長旺期，6 月中旬以前基本收割完畢，生長期大概 8 個月。本研究將 2020 年 10 月-2021 年 6 月(共 270 天)作為研究期，使用 DOW(Day of Wheat)表示研究期內的時間，表 1 為研究期內膠州市主要農作物的關鍵物候期。

　　1.3 衛星數據及預處理

　　Sentinel-2 是一種高分辨率多光譜成像衛星，攜帶有一枚多光譜成像儀(MSI)。它包括 Sentinel-2A 和 Sentinel-2B 兩顆衛星，一起提供 5d 間隔的影像。GEE 平臺提供了兩種級別的 Sentinel-2A/B 數據：Level-1C 和 Level-2A。Level-1C 級產品是指經過了輻射定標、幾何校正(主要包括正射校正、空間配準操作)的 TOA (Top-of-Atmosphere，大氣表層)反射率數據，Level-2A 級產品則是在 Level-1C 級產品的基礎上又經過了大氣校正的 SR(Surface Reflectance，地表反射率)數據[17]。本次研究使用 Sentinel-2 SR 數據集，已經消除了大氣吸收以及各項散射作用造成的誤差，能夠反演地物真實的表面反射率，使得作物分類更加準確。數據共含有 12 個光譜波段，分別具有不同的空間分辨率(10 m，20 m 和 60 m)，具體的數據介紹如表 2 所示。

　　研究期內共 218 景影像，由于原始數據是按 10 000 倍縮放的 SR 數據，將每景影像除以 10 000 來獲得真實單位的地表反射率數據。另外，數據還有三個 QA(Quality Assessment)波段，其中的 QA60 波段存儲有云掩膜信息，利用它消除了 Sentinel-2A/B 數據中被云和陰影污染的觀測，又使用線性插值的方法填補去云留下的缺失，以實現整個時間、空間域的完全覆蓋。

　　1.4 地面調查數據

　　樣本數據來源于 2021 年 3 月和 5 月膠州市野外實地調查數據，共獲得 436 個樣本點，其中冬小麥 212 個，非冬小麥 224 個。采樣時記錄樣本的經緯度及作物類別，并拍攝照片方便后續核查。

　　2 材料和方法

　　本研究主要包括以下兩個部分。首先利用圖像分割技術從高分辨率影像中提取膠州市的地塊，為了避免非農作物區域的干擾，將非耕地地塊剔除。然后在 GEE 平臺中使用實地調查得到的樣本結合構造的作物時序特征作為輸入，訓練基于 RF 算法的冬小麥分類器，將分類器分別應用于耕地地塊或像素，以實現膠州市冬小麥的提取。研究的技術路線如圖 1 所示。

　　2.1 耕地地塊獲取

　　在面向地塊的分類方法中地塊的獲取過程至關重要，地塊與真實地物的符合程度直接影響著信息提取的正確與否。圖像分割的過程其實是相鄰同質像素的結合和異質像素的分離過程，它在盡量減少圖像信息損失的基礎上將圖像分割成有意義的多邊形對象[18]。每個多邊形代表一個實地地物，內部具有相同的的光譜、紋理等信息。

　　本文使用易康(eCognition)軟件的多尺度分割方法來獲得耕地地塊，影像為下載得到的 2021 年 4 月的膠州市 GF2 PMS 數據，經過正射校正、波段融合后，空間分辨率達到了 1 m。為了使分割結果與影像更加吻合，經過多次試驗達到了最佳效果，這時分割尺度設置為 50，形狀因子為 0.3，緊致度為 0.4。采用目視解譯的方法將非耕地地塊(主要包括建筑物、道路等)刪除，只在耕地范圍內作圖，人工修正和細化邊界后得到了 99 733 個農作物地塊，局部效果如圖 2 所示。近年來國家耕地保護政策十分嚴格，獲得的耕地塊矢量具有較長期的穩定性，可以重復利用。

　　2.2 隨機森林分類

　　隨機森林(RF)算法是 Breiman 在 2001 年提出的一種強大的機器學習算法[19]。如今它已經在不同領域有了廣泛的應用，例如風險評估及預測、醫學檢驗和金融投資等[20-22]。它是基于決策樹弱分類器的集成學習算法，在利用集成算法優勢的同時規避了單個分類器的缺點，比許多傳統的分類器(比如單層神經網絡、最大似然法等)更具有準確性和魯棒性[23]。隨機森林的隨機主要有兩個方面：首先是對原始數據集采用bootstrap 方法(隨即且有放回的抽樣方法)抽取多個訓練集和其對應的測試集;其次是在生成決策樹時隨機抽取一部分特征用于樹節點的分裂。最后每個決策樹進行投票來確定最終分類結果。

　　RF 算法的模型可以表示為： ( ) arg Max I( ( ) ) n i i i H S h S Y ? ? ? ， (1) 式中， I 表示示性函數; Y 表示預測結果; ( ) i i h S 表示單顆決策樹。 RF 分類的步驟如下：(1) 從原始數據集 N 中隨機且有放回地抽取 n 個樣本，重復 n 次，而未抽中的樣本組成袋外數據集(OOB)，作為測試數據;(2) 從特征集合 M 中隨機抽取 m 個作為特征子集，且 m

　　2.2.1 光譜特征計算

　　NDVI(歸一化植被指數)是監測作物生長狀態的最佳指示因子，還能反映出植物冠層的背景(如土壤、潮濕地面、雪、枯葉等)影響，當土地還未耕作時裸露地表的 NDVI 指數值接近 0。NDVI 的計算公式為： NIR R NIR R NDVI ? ?? ????， (2) 式中，? NIR 為近紅外波段的反射率; R ?為紅波段的反射率。 EVI(增強植被指數)設置了更窄的紅邊波段，可以減少水汽的影響，加強了對植被稀疏地區的監測能力。另外它引入的藍光波段能夠降低土壤背景和大氣噪聲的干擾，可穩定地反映地表植被特征。EVI 的計算公式為： 2.5 6 7.5 1 R NIR R B NIR EVI ? ?? ? ??? ?? ? ? ? ?， (3) 式中， NIR ?為近紅外波段的反射率; R ?為紅波段的反射率， B ?為藍波段的反射率。

　　2.2.2 時序數據分析及合成

　　為了對比不同時間間隔的合成對時序曲線的影響，本文結合 Sentinel-2A/B NDVI、EVI 時間序列數據和實地調查的樣本數據，繪制不同作物原始的 NDVI、EVI 時間序列曲線。由于受到云、水汽、氣溶膠和傳感器等因素的干擾，時間序列數據容易出現異常值，時間序列曲線出現不規則波動現象，對波譜分析和特征構造造成不可忽略的影響[24]。平滑和去噪操作可以更加真實地反映各種作物的生長曲線，本文采用 S-G (Savitzky-Golay)濾波方法重構 Sentinel-2A/B 的時間序列數據，設置為移動窗口大小為 7(70 天觀測值)，多項式次數為 3。圖 3 為經過平滑后的作物 NDVI、EVI 時間序列曲線，橫坐標表示日期，單位(m/d)代表(年/月)，縱軸表示 NDVI、EVI 的值。

　　由 NDVI 時間序列曲線可以看出，冬小麥曲線整體呈現“兩峰一谷”的趨勢。夏玉米在 9 月下旬-10 月上旬收獲，然后冬小麥開始播種，7d 左右出苗，NDVI 逐漸增大;12 月份出現一個小的波峰，然后開始冬眠，光合作用衰減，NDVI 逐漸減小;波谷在 1 月份出現;在 2-3 月份開始返青起身，NDVI 迅速升高;4 月中下旬-5 月上旬進入抽穗期，達到第二個波峰;6 月成熟后收獲，植被指數值明顯下降。花生、春玉米和馬鈴薯都是春季作物，一般去年的夏季作物收獲之后都會為其保留耕地，NDVI 值較低，在 11 月中旬-次年 4 月都與冬小麥曲線有著明顯差異。EVI 曲線與 NDVI 較為類似，冬小麥曲線也有著“兩峰一谷”的形態，在 11 月中旬-次年 4 月都是區分冬小麥的關鍵時期。由于影像眾多，且直接使用單景影像易受到異常值的干擾，雖然進行了濾波處理但是作物曲線依然出現了振蕩現象。本文分別取 DOW：1-270 中的每 10 天、15 天和30 天的均值合成一景新的影像，就得到NDVI、 EVI 的不同時間間隔合成的時序影像。面向地塊的分類方法最后還要計算每個地塊內的均值作為地塊的特征。

　　2.2.3 特征選擇

　　經過以上計算后，得到了不同時間間隔合成的時序特征，如表 3 所示。

　　多光譜時序數據在提供更多光譜信息的同時，也帶來了高維輸入和輸出的新挑戰。由于大量光譜特征可能會攜帶高度相關的信息，增加模型的計算負擔，也許會出現過擬合現象。因此特征選擇就顯得十分重要，它在很大程度上決定著機器學習算法的效率。RF 允許使用基尼指數在分類時對特征重要性進行評價。因此使用 GEE 中的評價方法對特征波段進行排序，最后分別保留 14–24 個最佳特征進行冬小麥的提取。

　　2.2.4 模型訓練和分類

　　在進行地塊級別的分類時，由于膠州市面積較大，地塊數量較多，將它分為 3 個區分別計算和導出結果。面向地塊的分類方法將冬小麥樣本所屬的地塊標記為 1，非冬小麥地塊標記為 0;面向像素的分類方法將冬小麥樣本點標記為 1，非冬小麥點標記為 0。最后將樣本分別按照 7:3 的比例隨機分為訓練樣本和驗證樣本。將樣本和特征都導入 RF 分類器中進行分類，設置 RF 的參數如下：① numberOfTrees：樹數代表建立 RF 模型的決策樹的數量，在一定范圍內隨著樹數的增多，精度可能會略有提高，設置為 500;② seed：隨機種子設置為 999;③ 其他幾個參數：maxNodes(每棵樹中葉節點的最大數目)、minLeafPopulation(葉節點所需的最小樣本數)、bagFraction(每棵樹輸入到 bag 的分數)、variablesPerSplit(每次分割的變量數)保持默認設置。

　　為了方便與面向像素的分類方法進行對比，基于像素的冬小麥提取也使用多尺度分割得到的耕地地塊進行掩膜，時序特征及樣本也與面向地塊的分類方法保持一致。

　　2.3 結果評價指標

　　評價指標用來定量的衡量一個算法性能的優劣水平。混淆矩陣是最常用的方式，它能一目了然的展示出分類結果和地表真實信息的差異，如表 4 所示。

　　將真實類別為非小麥的樣本稱為負類，真實類別為小麥的樣本稱為正類，那么 TN(True Negative)代表負類的正確預測，FN(False Negative)代表正類的錯誤預測，FP(False Positive)代表負類的錯誤預測，TP(True Positive)代表正類的正確預測。樣本總和為 N，那么： N TN FN FP TP = ? ? ? 。 (4) 常用的幾種精度評價方式：總體分類精度(OA)、Kappa 系數、生產者精度(PA)和用戶精度(UA)都是基于混淆矩陣計算的。OA 表示正確分類的像元(地塊)總和除以總像元(地塊)數，計算公式如下： = TP TN TP TN OA TP TN FP FN N ? ??? ? ?。 (5) 但是 OA 只考慮分類正確的比例，不能對模型性能的好壞做出評價。Kappa 系數不僅考慮了被正確分類的像素(地塊)，又綜合了各種錯分和漏分誤差，是更為全面的的評價指標。Kappa 系數計算公式如下： 0 1 e e p p Kappa p ???， (6) 式中， 0 TP TN + p OA N ? ? ， (7) 2 ( )( ) ( )( ) e TN FP TN FN FN TP FP TP p N ? ? ? ? ?? 。 (8) 生產者精度是指分類器將所有像素(地塊)正確分為某類和某類真實總數的比例，例如冬小麥類別的生產者精度計算公式如下： TP PA TP FN ??。 (9) 用戶精度是指正確分到某類的像素(地塊)數與分類器分為某類的總數的比例，例如冬小麥類別的用戶精度計算公式如下：

　　3 結果與評價

　　3.1 時序曲線對比

　　圖 4 為冬小麥生長期內(DOW:0-270)不同時間間隔合成的作物生長曲線圖，橫軸表示時間，縱軸表示 NDVI、EVI 的值。

　　經過時間合成以后，曲線變得更加平滑也消除了一些異常波動，并且合成的時間間隔越長，作物曲線越平滑，但同時也失去了更多的細節，比如 30d 合成的冬小麥 NDVI、EVI 時序曲線均不再呈現“兩峰一谷” 的形態。時間間隔為 10 d 和 15d 的曲線既抵消了一些干擾也保留了作物生長的細節變化，在不同特征的曲線中冬小麥與其他作物的差異都比較明顯。

　　3.2 定性評價

　　為了檢驗本文方法的有效性，對比不同時序數據以及地塊和像素級別分類的優劣性，分別使用定性和定量的方式評價最終分類結果。定量方法即使用第 2 章中的評價指標進行對比分析，定性方法主要將提取結果與遙感影像進行對比，目視評判冬小麥提取效果，主要包括下面幾個方面：① 判斷冬小麥提取區域是否全面;② 判斷是否存在誤提現象;③ 基于像素的提取方法孤立像素現象是否嚴重。本文使用多尺度分割時利用的 GF2 1m 分辨率影像與各種方法得到的提取結果進行局部對比，圖 5 為典型區域的原始影像，冬小麥表現為深綠色，圖 6、圖 7、圖 8 分別為這些區域對應的提取結果。

　　在圖 6 中，30 d 和 15 d 合成的地塊、像素級別的提取均存在一定程度的誤提現象;10 d 合成的像素級別的提取有少部分的漏提;10 d 合成的地塊級別的提取結果符合影像顯示的冬小麥區域。在圖 7 中，地塊的不規整增大了提取的難度，基于像素的提取結果漏提現象嚴重;30 d 和 15 d 合成的地塊級提取也出現了漏提現象;10 d 合成的地塊級提取結果較為全面準確。圖 8 中，在復雜種植區域，基于像素的提取結果漏提和孤立像素現象都比較嚴重，提取效果不佳;30 d 合成的地塊級提取有小面積的漏提，15 d 和 10 d 合成的地塊級別的提取結果一致，在復雜種植區域的提取結果也較為理想。整體來看，地塊和像素級別的提取效果都會隨著合成時間序列數據時間間隔的縮小變得更好;當時間間隔相同時，地塊級別的提取并不全都優于像素級別的提取，但隨著時間間隔的縮小地塊級別的提取進步明顯，10 d 合成的面向地塊的提取方法總是最優的，幾乎沒有出現錯提、漏提現象。

　　3.3 定量評價

　　表 5 為最終得到的一系列量化結果，包括不同提取方法得到的 OA、Kappa 系數、冬小麥的 PA 和 UA。以上方法均得到了較高的精度，這證明了利用規律時間間隔合成的時間序列影像進行分類的有效性。整體而言，與基于像素的分類結果相比，面向地塊的分類方法都得到了更高的 OA 和 Kappa 系數，但 PA 和 UA 并不全都高于基于像素的分類。并且總體趨勢顯示，當合成 Sentinel-2A/B 影像的時間間隔越短時，冬小麥提取的精度越高。15 d 和 10 d 合成時間間隔相差不大，提取精度相差不大。

　　3.4 效率評價

　　表 6 為各種冬小麥提取方法所需要的計算時間和結果導出的時間。無論是計算時間還是導出時間，不同時間序列下地塊級別的提取之間或者像素級別的提取之間都相差不大，真正對提取時間產生影響的不是不同時間序列的合成方法，而是面向不同元素的分類方法。地塊與像素級別的提取在計算時間上最多相差 3 倍，但都在 1 分鐘之內，因此對計算效率的對比評價意義不大;然而在結果的導出上，兩者相差了 23 個小時，將近一天的時間，由于面向像素的分類方法精度也在比較高的水平，當在時間緊急的情況下，面向地塊的分類方法優勢不大。

　　4 冬小麥分布圖

　　通常在種植結構單一、農田連片區域的冬小麥更容易識別，但在地塊破碎區域由于內部的光譜變異和邊界的光譜混合導致識別較為困難。本文的研究區為整個膠州市，地塊數量眾多、類型多樣，方法具有普適性。由以上對比結果得知，10d 合成的面向地塊的分類方法精度更高，提取效果更好，因此利用此方法做出 2021 年的膠州市冬小麥種植分布圖，見圖 9，由圖可以明顯看出冬小麥種植區主要分布在膠州市的北部和西南部。

　　5 結束語

　　本文使用不同時間間隔合成的 Sentinel-2A/B 時間序列數據，提取了膠州市的冬小麥種植區，主要為了探索更適合提取冬小麥的時序數據與方法。為了對比在不同時序數據下面向像素和面向地塊的分類結果，兩種方法在完全相同的條件下進行。結果表明：① 當 Sentinel-2A/B 合成影像的時間間隔縮短時，可以為冬小麥與其它作物的區分提供更多信息，地塊和像素級別的提取準確度都會更高;② 在相同的時序數據下，面向地塊的分類方法并不總是優于面向地塊的分類方法，但它可以明顯改善面向像素分類結果中常見的“椒鹽現象”，且隨著合成時間序列數據時間間隔的縮小進步更明顯。同時，更短時間間隔(10 d)合成的 Sentinel-2 時序數據結合面向地塊的分類方法，對冬小麥的提取幫助更大，在一些不規整、復雜種植區域也獲得了很好的提取效果;③ 由于地塊級別的結果導出時間較長，當對精度和提取效果要求較高而時間比較充裕時，選擇地塊級別的分類更加合適;當對時間的要求較高而對精度的需求較低時，面向像素的分類更加符合條件。分類是不斷優化的過程，隨著遙感衛星研發技術、計算機和云平臺計算能力的不斷提高，現有的分類算法也在逐漸改進，未來更高級的圖像分割技術與更強大的云平臺計算能力的結合將獲得更加準確的冬小麥種植區提取結果。