基于Sentinel-2A遙感影像的濰坊市冬小麥種植面積提取研究*

孫逸飛，柳平增，張艷，宋長(zhǎng)青，張大雷，馬學(xué)文

(1. 山東農(nóng)業(yè)大學(xué)信息科學(xué)與工程學(xué)院，山東泰安，271000； 2. 山東農(nóng)業(yè)大學(xué)農(nóng)業(yè)大數(shù)據(jù)研究中心，山東泰安，271000；3. 農(nóng)業(yè)農(nóng)村部黃淮海智慧農(nóng)業(yè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，山東泰安，271000；4. 山東科技大學(xué)泰山科技學(xué)院，山東泰安，271000； 5. 山東農(nóng)大肥業(yè)科技有限公司，山東泰安，271000)

0 引言

冬小麥作為我國(guó)重要的糧食作物之一，精準(zhǔn)、高效的提取大范圍冬小麥種植面積和空間分布數(shù)據(jù)，可為農(nóng)業(yè)管理部門指導(dǎo)農(nóng)業(yè)生產(chǎn)、調(diào)整農(nóng)業(yè)結(jié)構(gòu)提供技術(shù)支撐[1]。濰坊市是我國(guó)重要的糧食生產(chǎn)基地，每年的小麥種植面積和產(chǎn)量都穩(wěn)居山東省前五位，近些年受設(shè)施農(nóng)業(yè)發(fā)展、城市化進(jìn)程加速等因素影響，對(duì)濰坊市冬小麥的種植產(chǎn)生了較大的沖擊。因此，精準(zhǔn)、高效地獲取冬小麥種植面積、空間分布以及變化趨勢(shì)，對(duì)于提高冬小麥針對(duì)性管理的高效性、種植結(jié)構(gòu)布局的有效性等有重要的推動(dòng)作用[2]。但傳統(tǒng)的基于抽樣和統(tǒng)計(jì)部門逐級(jí)上報(bào)農(nóng)情的調(diào)查方法，效率低、結(jié)果不直觀且容易出現(xiàn)人為因素造成的誤差[3]。遙感技術(shù)具有覆蓋范圍廣、重訪周期短，獲取成本相對(duì)低等優(yōu)勢(shì)，在作物長(zhǎng)勢(shì)監(jiān)測(cè)、農(nóng)業(yè)氣象以及產(chǎn)量估算中扮演著重要的角色，且對(duì)大面積露天農(nóng)業(yè)生產(chǎn)的調(diào)查、評(píng)價(jià)、監(jiān)測(cè)和管理具有獨(dú)特的作用[4-5]。

隨著遙感技術(shù)的發(fā)展，中高空間分辨率的遙感衛(wèi)星影像數(shù)據(jù)資源增多，并開始逐漸應(yīng)用于農(nóng)作物空間分布信息提取。但中高空間分辨率遙感影像細(xì)節(jié)更加豐富，傳統(tǒng)的基于單像元的分類方法，已經(jīng)難以滿足實(shí)際應(yīng)用的精度要求。近年來(lái)，國(guó)內(nèi)外學(xué)者針對(duì)面向?qū)ο?、深度學(xué)習(xí)方法在高效精確的目標(biāo)地物提取上進(jìn)行了積極探索。金永濤等[6]利用深度學(xué)習(xí)結(jié)合面向?qū)ο髮?duì)典型地物進(jìn)行提取，提取精度較為理想。白宇[7]利用深度學(xué)習(xí)算法識(shí)別不同地區(qū)的林地范圍。陳燕生等[8]使用基于U-Net自主改進(jìn)的卷積神經(jīng)網(wǎng)絡(luò)，對(duì)小宗作物遙感圖像進(jìn)行高精度分割研究，實(shí)現(xiàn)精細(xì)化分割效果。Rezaee等[9]使用CNN分類框架對(duì)不同類型的濕地進(jìn)行識(shí)別，試驗(yàn)結(jié)果表明CNN在復(fù)雜濕地分類提取中效果優(yōu)于隨機(jī)森林。Duarte-Carvajalino等[10]采用淺層學(xué)習(xí)和深度學(xué)習(xí)方法對(duì)土豆作物的生長(zhǎng)狀況進(jìn)行監(jiān)測(cè)識(shí)別，試驗(yàn)結(jié)果表明CNN和RF的效果優(yōu)于DNN、SVM。

本研究利用Sentinel-2A遙感影像數(shù)據(jù)，以濰坊市為研究區(qū)，在實(shí)地調(diào)研掌握研究區(qū)域地物分布真實(shí)樣本基礎(chǔ)上，利用面向?qū)ο笈c深度學(xué)習(xí)相結(jié)合的方法，綜合利用地物光譜及紋理等特征，選取合適的分割尺度，訓(xùn)練出最優(yōu)冬小麥提取模型，利用混淆矩陣對(duì)分類結(jié)果進(jìn)行精度評(píng)價(jià)。

1 研究區(qū)概況、數(shù)據(jù)源、預(yù)處理及技術(shù)路線

1.1 研究區(qū)概況

濰坊市(118°10′E～120°01′E，35°41′N～37°26′N)位于山東半島中部，南依泰沂山脈，北瀕渤海萊州灣，市域自北向南，由低到高，大體上分為低地、平原及低山丘三個(gè)地貌區(qū)。濰坊市域?qū)儆诒睖貛Ъ撅L(fēng)區(qū)，受典型季風(fēng)的季風(fēng)氣候的影響，四季分明，氣候溫和，雨量適中，光照充足，適宜農(nóng)作物的種植。冬春季種植作物主要有冬小麥、玉米、棉花、油菜等，冬小麥和玉米是該地主要的糧食作物。濰坊市作物結(jié)構(gòu)較為復(fù)雜，地塊破碎，使得精確提取冬小麥種植空間分布具有一定的挑戰(zhàn)性?？焖偬崛《←湻N植區(qū)域信息，分析2017—2021年冬小麥種植的空間分布特征，對(duì)濰坊市農(nóng)業(yè)產(chǎn)業(yè)政策的制定及結(jié)構(gòu)調(diào)整起到積極的參考作用。

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

1.2.1 Sentinel-2A數(shù)據(jù)獲取

本研究綜合考慮影像質(zhì)量、作物生長(zhǎng)和分布情況等因素，選擇Sentinel-2A衛(wèi)星影像數(shù)據(jù)。Sentinel-2A衛(wèi)星設(shè)計(jì)了一種可靠的多光譜地球觀測(cè)系統(tǒng)，將高空間分辨率與寬視場(chǎng)結(jié)合，并采用三軸姿態(tài)控制，能夠獲得優(yōu)越的精度和穩(wěn)定性。重訪周期為10 d，并攜帶覆蓋13個(gè)波段的多光譜成像儀，除可見光及近紅外、短波紅外等常用波段外，還包括3個(gè)植被生長(zhǎng)與變化敏感的紅邊波段(表1)，紅邊波段的引入增加了地物之間的可分性測(cè)度，減少了影像噪聲造成的景觀破碎度的增加，有利于作物種植面積提取[11]。

表1 Sentinel-2A影像多光譜波段主要參數(shù)信息

在中華人民共和國(guó)農(nóng)業(yè)農(nóng)村部種植業(yè)管理司獲取到濰坊市主要作物物候期表(表2)，由表2可知從種植當(dāng)年12月到次年3月是冬小麥的苗期，而其他大宗農(nóng)產(chǎn)品尚未播種或返青，是識(shí)別冬小麥的關(guān)鍵生育期。通過(guò)歐空局官方網(wǎng)站(https://scihub.copernicus.eu/)獲取2017年—2021年2月到3月份L1C級(jí)別的濰坊市Sentinel-2A遙感數(shù)據(jù)(影像云量小于5%)。

表2 濰坊市主要作物物候期表

1.2.2 其他資料

1) 濰坊市各縣區(qū)行政區(qū)劃矢量數(shù)據(jù)，數(shù)據(jù)來(lái)源于全國(guó)地理信息資源目錄服務(wù)系統(tǒng)。

2) 濰坊市2017—2021年統(tǒng)計(jì)年鑒和實(shí)際調(diào)研數(shù)據(jù)。

1.2.3 Sentinel-2A數(shù)據(jù)預(yù)處理

影像在成像過(guò)程中會(huì)受到衛(wèi)星速度變化、電磁場(chǎng)與大氣相互作用、隨機(jī)噪聲等因素的干擾，地物波譜信息發(fā)生變化，為改善影像的解譯條件，提高可辨性，需對(duì)影像進(jìn)行預(yù)處理。預(yù)處理包括：輻射定標(biāo)、大氣校正、像素重采樣、鑲嵌裁剪、波段融合等操作。

1.3 技術(shù)路線

本研究基于Sentinel-2A影像，對(duì)影像進(jìn)行預(yù)處理后，分別利用面向?qū)ο蠓椒ê碗S機(jī)森林方法進(jìn)行試驗(yàn)并對(duì)比精度。選擇分類精度較高的算法進(jìn)一步進(jìn)行深度學(xué)習(xí)模型訓(xùn)練。具體流程見圖1。

圖1 總體流程

2 基于面向?qū)ο笈c隨機(jī)森林分類方法提取

2.1 面向?qū)ο蠓椒?/h3>

面向?qū)ο筇崛×鞒倘鐖D2所示。

圖2 面向?qū)ο筇崛×鞒?/p>

面向?qū)ο蟮倪b感分類方法是對(duì)利用遙感數(shù)據(jù)的空間、紋理和光譜特征建立判別規(guī)則或判別函數(shù)的推理對(duì)分割后的目標(biāo)對(duì)象進(jìn)行分類提取[12-13]。該方法能有效避免“同物異譜”和“同譜異物”的問(wèn)題，避免出現(xiàn)“椒鹽現(xiàn)象”[14-15]。關(guān)鍵技術(shù)包括：影像分割、影像合并、特征選擇和規(guī)則建立。

2.1.1 影像分割與影像合并

采用多尺度分割算法(Multiresolution Segmentation)實(shí)現(xiàn)影像的分割。多尺度面向?qū)ο蠓指钏惴▽D像看作是一個(gè)由區(qū)域和區(qū)域之間的拓?fù)潢P(guān)系組成的一張區(qū)域臨接圖[16]，根據(jù)指定尺度進(jìn)行分割，采用從單像元大小的區(qū)域開始相鄰影像區(qū)域兩兩合并增長(zhǎng)的方法，設(shè)定閾值控制合并區(qū)域，保證生成高度同質(zhì)性(或異質(zhì)性最小)的影像分割區(qū)域(影像對(duì)象)[17]，從而適應(yīng)于最佳分離和表示地物目標(biāo)。在進(jìn)行影像區(qū)域的合并操作時(shí)，以區(qū)域鄰接圖中的各個(gè)區(qū)域?yàn)檠芯繉?duì)象，找到滿足局部最優(yōu)合并條件的區(qū)域?qū)?，將最?yōu)合并條件區(qū)域?qū)M(jìn)行合并操作，與此同時(shí)更新區(qū)域?qū)υ鞠噙B接的區(qū)域的特征值，及其與更新區(qū)域的合并權(quán)重。

為了得到適合冬小麥提取的最佳分割參數(shù)，基于Sobel算子的邊緣影像分割尺度、Full Lambda-Schedule算法影像片段合并和紋理內(nèi)核參數(shù)進(jìn)行了大量試驗(yàn)。最終選取冬小麥影像錯(cuò)分率最小的分割尺度。影像最終分割效果圖如圖3所示。

圖3 分割效果圖

2.1.2 特征選擇、分類規(guī)則建立

影像分割完成后，多邊形對(duì)象包括很多屬性信息。主要包括：光譜特征、形狀特征、紋理特征、顏色空間和波段比值等[18]。將分割形成的有特定含義的子區(qū)域定義為一個(gè)個(gè)的影像對(duì)象，每個(gè)對(duì)象都有其特有的屬性信息，存放在矢量面文件中。將屬性信息進(jìn)行進(jìn)一步挖掘，提取出特征變量。采用基于規(guī)則特征的面向?qū)ο蟮奶崛》椒?，利用目?biāo)地物的波譜、紋理、RENDVI值、形狀、顏色的差異，進(jìn)行冬小麥區(qū)域的提取。圖4為組合分類特征后冬小麥提取效果圖細(xì)節(jié)展示。

圖4 冬小麥種植面積提取效果圖(面向?qū)ο?

2.2 隨機(jī)森林分類

隨機(jī)森林技術(shù)路線圖如圖5所示。

圖5 隨機(jī)森林算法技術(shù)路線圖

隨機(jī)森林通過(guò)對(duì)數(shù)據(jù)及特征變量進(jìn)行隨機(jī)重采樣，構(gòu)建多個(gè)CART決策樹，最終采用多決策樹投票的方式確定數(shù)據(jù)的類別歸屬。它能夠處理具有高維特征的輸入樣本，又對(duì)過(guò)度擬合不敏感[19]。在生成過(guò)程中可以對(duì)誤差建立一個(gè)無(wú)偏估計(jì)，因此對(duì)于遙感影像農(nóng)作物面積提取有良好的抗噪性能，在農(nóng)業(yè)制圖研究方面取得了良好的分類結(jié)果[20-21]。隨機(jī)森林算法可以處理具有多維特征的輸入樣本，并且不需要降維，在生成過(guò)程中能夠獲取到內(nèi)部生成誤差的一種無(wú)偏估計(jì)，對(duì)于缺省問(wèn)題也能獲得良好的結(jié)果，不需要做很多參數(shù)的調(diào)試。設(shè)置隨機(jī)森林分類器參數(shù)max_features取值為6，選擇Gini Coefficient函數(shù)。number_of_trees取值為100。

圖6為本次分類中選擇的訓(xùn)練樣本樣例，經(jīng)分類后，輸出此類分類結(jié)果，結(jié)果如圖7所示。

(a) 冬小麥

圖7 冬小麥提取效果圖(隨機(jī)森林)

2.3 分類結(jié)果精度對(duì)比

隨機(jī)選取研究影像中的500個(gè)樣本點(diǎn)像元，通過(guò)混淆矩陣計(jì)算出總體精度和Kappa系數(shù)，對(duì)分類結(jié)果進(jìn)行精度評(píng)價(jià)，結(jié)果如圖8所示。結(jié)果表明面向?qū)ο蟮姆诸惙椒ńY(jié)果優(yōu)于隨機(jī)森林算法。故接下來(lái)的研究采用面向?qū)ο蠓椒ǖ贸龅姆诸惤Y(jié)果創(chuàng)建標(biāo)簽進(jìn)行深度學(xué)習(xí)。

圖8 精度評(píng)價(jià)

3 基于深度學(xué)習(xí)算法冬小麥種植面積提取模型

3.1 深度學(xué)習(xí)算法與U型網(wǎng)絡(luò)

深度學(xué)習(xí)是機(jī)器學(xué)習(xí)的分支，它可以通過(guò)學(xué)習(xí)一種深層非線性網(wǎng)絡(luò)結(jié)構(gòu)，實(shí)現(xiàn)復(fù)雜函數(shù)的逼近，展現(xiàn)了強(qiáng)大的從少數(shù)樣本中集中學(xué)習(xí)數(shù)據(jù)集本質(zhì)特征的能力[22]。ENVINet5模型基于Tensorflow框架開發(fā)，模型架構(gòu)基于U-Net神經(jīng)網(wǎng)絡(luò)改進(jìn)。該架構(gòu)基于Ronneberger、Fischer和Brox(2015)開發(fā)的U-Net架構(gòu)，用于對(duì)圖像中的每個(gè)像素進(jìn)行分類[23]。

圖9為ENVINet5模型架構(gòu)。

圖9 ENVINet5模型架構(gòu)

該架構(gòu)有5個(gè)級(jí)別和27個(gè)卷積層，每個(gè)級(jí)別表示模型中不同的像素分辨率。架構(gòu)的上下文視場(chǎng)表示在訓(xùn)練期間周圍區(qū)域?qū)γ總€(gè)像素的貢獻(xiàn)程度。將原始遙感圖像輸入神經(jīng)網(wǎng)絡(luò)模型，進(jìn)行切片，切片大小大于上下文視場(chǎng)時(shí)允許一次執(zhí)行更多的訓(xùn)練，并更快的分類。再通過(guò)3×3的卷積層對(duì)切片進(jìn)行卷積，使得圖像特征增多，使用Rule激活函數(shù)避免訓(xùn)練中出現(xiàn)過(guò)擬合現(xiàn)象。一個(gè)epoch表示在偏倚判定調(diào)整之前訓(xùn)練的切片數(shù)。通過(guò)2×2池化層減小圖像尺寸，保留重要特征。使用逆變換抽樣的統(tǒng)計(jì)技術(shù)，使得傳入模型的樣本與它們對(duì)概率密度函數(shù)的貢獻(xiàn)成正比。加入偏差損失函數(shù)會(huì)使得損失函數(shù)偏向于尋找特征像素。提供一個(gè)實(shí)體距離Solid Distance參數(shù)來(lái)擴(kuò)展線狀特征大小，使得能夠完整地表示冬小麥區(qū)域。為解決深度學(xué)習(xí)難以學(xué)習(xí)銳利邊緣的問(wèn)題，引入Blur Distance參數(shù)，設(shè)置一個(gè)模糊邊緣。

3.2 提取過(guò)程

TensorFlow模型是由一組基本的神經(jīng)網(wǎng)絡(luò)參數(shù)定義的，為了讓模型能夠提取特定目標(biāo)，必須先使用一組能夠指明這些特征的標(biāo)簽圖像對(duì)模型進(jìn)行訓(xùn)練，提取結(jié)果為類激活灰度圖CAM，在類激活灰度圖中，DN值代表該像元屬于目標(biāo)的概率。該方法在ENVI5.5中使用Deeplearning模塊實(shí)現(xiàn)，操作步驟主要有：創(chuàng)建標(biāo)簽柵格、訓(xùn)練模型、模型分類、得到分類結(jié)果。技術(shù)路線如圖10所示。

圖10 基于深度學(xué)習(xí)算法冬小麥種植面積提取流程

3.2.1 創(chuàng)建標(biāo)簽圖像

在創(chuàng)建標(biāo)簽圖像時(shí)，需要輸入一個(gè)收集感興趣區(qū)域的樣本。目標(biāo)特征分為一維、二維和點(diǎn)狀。支持點(diǎn)、像元、多邊形、折線或閾值。采用二維特征標(biāo)記冬小麥的輪廓。標(biāo)簽圖像包含輸入圖像的原始波段和一個(gè)掩膜波段。掩膜波段DN值為1的像元代表目標(biāo)，0代表背景。

3.2.2 初始化及訓(xùn)練TensorFlow模型

在初始化TensorFlow模型時(shí)，需定義該模型的體系結(jié)構(gòu)、切片大小(patch size)、訓(xùn)練所用波段等。訓(xùn)練是將標(biāo)簽圖像反復(fù)暴露在模型中。隨著訓(xùn)練的進(jìn)行，模型將學(xué)習(xí)把標(biāo)簽圖像中的光譜和空間信息轉(zhuǎn)換為CAM灰度圖。在訓(xùn)練過(guò)程中，模型嘗試初始猜測(cè)并生成一個(gè)隨機(jī)CAM灰度圖，將其與標(biāo)簽圖像中的掩膜波段進(jìn)行對(duì)比。通過(guò)擬合優(yōu)度函數(shù)(損失函數(shù))可以識(shí)別隨機(jī)猜測(cè)的結(jié)果哪些是錯(cuò)誤的。通過(guò)調(diào)整模型的權(quán)重，對(duì)提取效果進(jìn)行優(yōu)化。

3.2.3 執(zhí)行分類及批處理操作

將訓(xùn)練好的模型在其他圖像中定位相同的特征。并將2017—2021年研究區(qū)域內(nèi)影像在該模型中進(jìn)行批處理分類操作，圖11是深度學(xué)習(xí)冬小麥種植面積提取結(jié)果圖。

圖11 基于深度學(xué)習(xí)算法冬小麥種植面積提取效果圖

3.2.4 分類后處理

將相同年份的影像進(jìn)行鑲嵌、裁剪、小圖斑處理，并進(jìn)行統(tǒng)計(jì)，得出2017—2021年濰坊市冬小麥種植面積并與統(tǒng)計(jì)年鑒中抽樣調(diào)查得到的濰坊市冬小麥種植面積數(shù)據(jù)進(jìn)行對(duì)比，如表3所示。

表3 濰坊市冬小麥種植面積提取結(jié)果

從2017—2021年濰坊市冬小麥種植面積提取結(jié)果可以看出，冬小麥主要分布在濰坊市的東南部，諸城市和高密市為濰坊市冬小麥種植面積較大的兩個(gè)地級(jí)市，高密市的冬小麥種植密度最高。2017—2021年間的冬小麥種植面積變化于降低的狀態(tài)，面積降低約7.4%。

3.3 精度評(píng)價(jià)

利用高分辨率的Google影像建立真實(shí)的冬小麥樣本區(qū)域，對(duì)濰坊市2017—2021年冬小麥提取結(jié)果分別進(jìn)行混淆矩陣精度評(píng)價(jià)并取平均。并加入制圖精度、用戶精度兩個(gè)評(píng)價(jià)指標(biāo)。

制圖精度(Producer Accuracy)表示實(shí)際中某個(gè)地物i樣本在分類結(jié)果中被正確分類為地物i的概率。制圖精度計(jì)算公式

(1)

用戶精度(User Accuracy)表示從分類結(jié)果中任意取一個(gè)地物i像元，其在實(shí)際中是地物i類型的概率。用戶精度計(jì)算公式

(2)

式中：Pii——混淆矩陣第i行第i列的像元數(shù)；

P+i、Pi+——混淆矩陣第i行和第i列的像元總數(shù)。

2017—2021年深度學(xué)習(xí)冬小麥提取制圖精度達(dá)到97.13%，用戶精度達(dá)到97.87%，總體精度達(dá)到93.1%，誤差精度2.81%，Kappa系數(shù)0.91，該精度滿足以市級(jí)或直轄市為主的局部尺度的農(nóng)情遙感監(jiān)測(cè)精度的需求。

4 結(jié)論

本研究在綜合考慮冬小麥大尺度種植面積監(jiān)測(cè)業(yè)務(wù)運(yùn)行成本和信息提取精度的前提下，利用圖像分割、RENDVI、紋理特征作為分類特征，分別利用面向?qū)ο蟆㈦S機(jī)森林方法進(jìn)行提取，并利用精度較高的方法結(jié)果建立標(biāo)簽圖像，再利用深度學(xué)習(xí)技術(shù)建立了冬小麥提取模型，分別提取了2017—2021年濰坊市冬小麥種植面積。利用ENVI和ArcGIS統(tǒng)計(jì)各縣區(qū)冬小麥種植提取面積，并與統(tǒng)計(jì)調(diào)查數(shù)據(jù)進(jìn)行精度檢驗(yàn)，最終獲得冬小麥提取結(jié)果。得出以下結(jié)論。

1) 通過(guò)RENDVI時(shí)間序列及反射光譜曲線，選擇冬小麥遙感提取的最佳時(shí)相，相比于NDVI時(shí)間序列可以更好地將大蒜與冬小麥區(qū)分開來(lái)。利用圖像分割技術(shù)并加入紋理特征后，面向?qū)ο蠓诸惖奶崛【认啾扔陔S機(jī)森林方法有顯著提高。

2) 在利用深度學(xué)習(xí)建立冬小麥提取模型時(shí)，對(duì)圖像進(jìn)行隨機(jī)裁剪和全景縮放，歸一化以及閾值處理并利用監(jiān)督分類結(jié)合手工采集來(lái)解決標(biāo)簽冬小麥數(shù)據(jù)集采集分類問(wèn)題，提高了模型的訓(xùn)練速度和精度。與面向?qū)ο蠓诸惙椒ㄏ啾?，深度學(xué)習(xí)方法精度更高，總體精度提高了4%，Kappa系數(shù)提高了0.07，且有效的減少了背景地物的干擾。