哨兵2號多時(shí)相植被指數(shù)作物分類及監(jiān)測

祁向前,孫德浩,賈連星

(1.龍巖學(xué)院 資源工程學(xué)院，福建 龍巖 364012;2.中國礦業(yè)大學(xué) 環(huán)境與測繪學(xué)院,江蘇 徐州 221116;3. 黑龍江科技大學(xué) 礦業(yè)工程學(xué)院,哈爾濱 150022)

以往對農(nóng)作物的監(jiān)測(分類、面積變化等)大多靠人工現(xiàn)場采集信息，既耗時(shí)又費(fèi)力，近年來隨著現(xiàn)代化腳步的加快，多種類的衛(wèi)星可以滿足不同的監(jiān)測需求，分辨率的提高也保障了監(jiān)測信息的準(zhǔn)確性，衛(wèi)星用于監(jiān)測農(nóng)業(yè)生產(chǎn)，能夠體現(xiàn)其快速、客觀、準(zhǔn)確的特點(diǎn)，有利于及時(shí)掌握農(nóng)情基本信息，進(jìn)行宏觀調(diào)控。國內(nèi)外學(xué)者研究表明，利用多時(shí)相遙感影像植被指數(shù)和不同作物的物候特征結(jié)合進(jìn)行作物識別與監(jiān)測是一種行之有效的方法。例如Jakubauskas等[1]采用諧波算法對NOAA/AVHRR-NDVI時(shí)序數(shù)據(jù)進(jìn)行研究，得到作物分類信息。Moody[2]等利用歸一化差值植被指數(shù)(Normailized Difference Vegetation Index,NDVI )時(shí)間序列的Fourier變換來劃分土地的覆蓋類型等。Lanjeri[3]等利用多時(shí)相掩膜法成功提取了葡萄園的種植面積。競霞等[4]利用多時(shí)相Landsat TM數(shù)據(jù)，提取冬小麥不同時(shí)相的作物波譜信息，成功監(jiān)測了北京地區(qū)的冬小麥種植結(jié)構(gòu)調(diào)整情況。Young和[5]Running[6]通過分析研究區(qū)植被類型的年NDVI變化情況，并結(jié)合作物的物候特征對植被類型進(jìn)行分類。Conese等[7]對遙感數(shù)據(jù)進(jìn)行了單時(shí)相和多時(shí)相的研究，發(fā)現(xiàn)多時(shí)相數(shù)據(jù)能夠極大改善分類精度。賈樹海等[8]針對3個(gè)不同時(shí)期的遙感影像，結(jié)合不同作物的物候特征提取其NDVI植被指數(shù)，采用監(jiān)督分類的方法，較好地提取了花生的種植面積并制圖。

哨兵2號是由兩顆極軌衛(wèi)星組成的星座，具有高分辨率多光譜成像的特點(diǎn)，一顆衛(wèi)星的重訪周期為10 d，兩顆互補(bǔ)，重訪周期為5 d，在紅邊范圍內(nèi)有較多波段，有利于對植被的監(jiān)測。國內(nèi)外均有學(xué)者利用單時(shí)相哨兵2號數(shù)據(jù)實(shí)現(xiàn)了作物識別[9-10]。NDVI作為一種植被指數(shù)，被廣泛應(yīng)用于反映作物在不同季節(jié)、不同生育期表現(xiàn)出來的生理特征，也有效的應(yīng)用于土地利用覆蓋監(jiān)測、植被覆蓋密度評價(jià)、作物識別和作物產(chǎn)量預(yù)報(bào)等[11-13]。哨兵數(shù)據(jù)的波段最高分辨率為10 m，單純依靠目視解譯對作物進(jìn)行分類較為困難，因此構(gòu)建多時(shí)相NDVI植被指數(shù)來研究就能夠很好地提高分類精度[14-16]。

上述文獻(xiàn)主要是針對不同種類農(nóng)作物進(jìn)行分類或者是對單一農(nóng)作物進(jìn)行面積提取，涉及到的農(nóng)作物精細(xì)分類(例如夏玉米、春玉米、山地夏玉米等)及多種農(nóng)作物種植面積變化研究還很少。文中以河北省石家莊市靈壽縣的多種作物及地貌類型為例，采用多時(shí)相決策樹分類法，利用所建立的作物及地貌樣本點(diǎn)不同時(shí)期內(nèi)監(jiān)測到的NDVI值建立時(shí)序NDVI曲線，并在具有明顯區(qū)分特征的月份確定閾值，對研究區(qū)域內(nèi)的多種作物及地貌類型進(jìn)行分類研究。

1 研究區(qū)與數(shù)據(jù)處理

1.1 研究區(qū)

靈壽縣(113°45′E～114°28′E，38°16′N～38°48′N)位于河北省西部，太行山東側(cè)，在省會石家莊西北35 km處，東南與正定縣毗連，南與鹿泉區(qū)隔滹沱河相望，西同平山縣接壤，西北隅與平山縣、阜平縣及山西省五臺縣接境(圖1)。全境呈西北—東南走向，寬15～20 km，長約100 km，靈壽縣屬北溫帶亞濕潤氣候，處于半干旱、半濕潤季風(fēng)區(qū)。年平均氣溫13.7 ℃，無霜期190 d左右，年平均降雨量417 mm。全縣總面積1 066.2 km2，耕地面積521.12 km2，約占全縣總面積的48.87%，主要分布在中南部區(qū)域，雙熟制耕地主要有冬小麥-花生和冬小麥-玉米大豆，單熟制耕地主要有春玉米、夏玉米、山地夏玉米，縣域的北邊多為山區(qū)林地，其他的植被及地貌類型有雜草、河流湖泊和居民地道路裸地等。土壤有黏土、沙壤等多種類型土質(zhì)，河流主要有滹沱河、慈河、松陽河、衛(wèi)水河等。

圖1 2020年研究區(qū)域

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

根據(jù)作物的物候特征了解到作物有明顯植被指數(shù)區(qū)分度的月份為4月、6月、8月和10月，后續(xù)用所監(jiān)測樣本點(diǎn)的時(shí)序植被指數(shù)變化分析也驗(yàn)證了該想法，并且考慮到云量會對影像分類精度造成影響，因此通過綜合比較各期影像云量覆蓋率，分別選取2019—2020年兩年的8期影像，每年4期，分別是4月、6月、8月和10月(表1)，涵蓋了全部農(nóng)作物的主要生長季(其中4月、6月、10月影像無云量覆蓋，8月影像云量覆蓋度為0.8%)。8期影像數(shù)據(jù)均來自歐洲航天局(ESA)(https://scihub.copernicus.eu/dhus/#/home)，哨兵2號由哨兵2A和哨兵2B兩顆衛(wèi)星組成，重訪周期為5 d，共由13個(gè)波段組成，其中較多的紅外波段有利于植被的監(jiān)測，分辨率有10 m、20 m、60 m，經(jīng)過處理后可達(dá)到10 m。網(wǎng)站的影像數(shù)據(jù)出現(xiàn)離線狀態(tài)時(shí)將其收藏并發(fā)送下載請求，過半個(gè)至一個(gè)工作日即可下載。

表1 影像清單及各農(nóng)作物發(fā)育期

從歐空局下載下來的數(shù)據(jù)是L1C級，無需進(jìn)行輻射定標(biāo)，經(jīng)過大氣校正后進(jìn)行重采樣(10 m)，進(jìn)行波段融合，剔除掉不需要的海岸/氣溶膠波段(band.1)和水蒸氣波段(band.9)。經(jīng)過裁剪和空間配準(zhǔn)，得到處理后的影像數(shù)據(jù)。

圖2 技術(shù)流程圖

2 方法研究

2.1 技術(shù)路線

研究技術(shù)路線如圖2所示，在2019年3月通過外業(yè)GPS定位點(diǎn)來建立樣本點(diǎn)監(jiān)測NDVI值的變化(圖3)，從北至南分別為雜草、林地、山地-夏玉米、河流、湖泊、春玉米、夏玉米、冬小麥-花生、居民地、道路、裸地，冬小麥-玉米大豆，監(jiān)測時(shí)間持續(xù)到10月底所有作物收割結(jié)束(或正在收割)的NDVI值的變化情況，繪制NDVI時(shí)間序列曲線。了解農(nóng)作物的物候特征，選取具有代表性的8期多時(shí)相遙感數(shù)據(jù)，預(yù)處理后計(jì)算每幅影像的歸一化植被指數(shù)NDVI，NDVI算式如式(1)。根據(jù)農(nóng)作物多時(shí)相NDVI變化情況建立分類判別規(guī)則，確定相應(yīng)的閾值將各種作物及其他植被和地貌進(jìn)行分類并進(jìn)行分類統(tǒng)計(jì)，得到兩年以來各種作物的種植面積的變化情況。為了驗(yàn)證分類結(jié)果的準(zhǔn)確性，于2019年8月及2020年8月到現(xiàn)場勘察，采集多個(gè)樣本點(diǎn)建立驗(yàn)證樣本集，對當(dāng)年的作物及其他地貌分類通過混淆矩陣進(jìn)行精度評價(jià)來驗(yàn)證分類的準(zhǔn)確性。

(1)

式中:BNIR為哨兵2號近紅外波段(第八波段)的地表反射率;BRED為哨兵2號紅光波段(第四波段)的地表反射率。

圖3 樣本點(diǎn)位置

2.2 植被指數(shù)曲線及可分性

本研究區(qū)域內(nèi)典型的地表類型有雙熟制耕地(冬小麥-花生、冬小麥-玉米大豆)，單熟制耕地(春玉米、夏玉米、山地-夏玉米)，林草類(林地、雜草)，非植被地表(居民地道路裸地、河流湖泊)。其中非植被地表的居民地道路裸地的植被指數(shù)常年較低，在0～0.1附近，河流湖泊的植被指數(shù)全年低于0，在-0.1～0之間。雙熟制耕地NDVI曲線具有雙峰，根據(jù)4月影像可以與單熟制作物區(qū)分開，10月底具有較低植被指數(shù)可與林地做區(qū)分，花生收割時(shí)間相比大豆較早，10月底花生收割完(植被指數(shù)低于0.20)，大豆正在收割(植被指數(shù)大于0.20)。單熟制耕地中，春玉米4月份播種，6月與其它兩種玉米相比有較高的植被指數(shù)，夏玉米和山地-夏玉米6月播種，此時(shí)植被指數(shù)很低，三者在8月植被指數(shù)到達(dá)峰值，其中山地-夏玉米由于多種植在山地丘陵地區(qū)，自然條件惡劣且密度低，植被指數(shù)不超過0.6。林草類植被由于不是人工種植，在8月后植被指數(shù)并不會出現(xiàn)大降幅，其中林地常年較高，雜草多出現(xiàn)在道路兩旁及山區(qū)，植被指數(shù)在8月到達(dá)峰值時(shí)低于作物的峰值(圖4)。

圖4 農(nóng)作物及林草類NDVI時(shí)序曲線

在構(gòu)建的NDVI時(shí)序曲線中，非植被地表由于全年低植被指數(shù)因此有著極高的識別度；雙熟制耕地作物及林地特征差異明顯，識別度較高；3種玉米具有相同的植被指數(shù)變化趨勢，在起伏及峰值上有差異，也較容易識別；雜草種類多樣，密度低，分布不集中，不能獲得全部種類雜草的植被指數(shù)變化情況，只能得到一個(gè)大致變化范圍，識別度一般。綜上所述，利用NDVI植被指數(shù)來對作物及其他地表類型進(jìn)行識別分類具有較好的可行性。

2.3 決策樹方法用于作物及地貌識別

決策樹作為一種分類方法，具有靈活、直觀、運(yùn)算效率高的優(yōu)點(diǎn)。它是一種從無次序、無規(guī)則的樣本數(shù)據(jù)集中推理出決策樹表示形式的分類規(guī)則方法，采用自頂向下的遞歸方式，在決策樹的內(nèi)部節(jié)點(diǎn)進(jìn)行屬性值的比較并根據(jù)不同的屬性值判斷從該節(jié)點(diǎn)向下的分支，在決策樹的葉節(jié)點(diǎn)得到結(jié)論。該方法在遙感影像分類和信息提取中已得到廣泛應(yīng)用。文中根據(jù)上述分析，設(shè)定不同的閾值來構(gòu)建決策樹對作物及地貌進(jìn)行分類研究。

閾值的設(shè)置必須滿足具有明顯區(qū)分度的要求，才能保證較高精度的分類。分別將4月、6月、8月、10月影像的NDVI值設(shè)為B1、B2、B3、B4。根據(jù)參考時(shí)序曲線以及多次改進(jìn)閾值可以得到，當(dāng)B2大于0.6且B4大于0.5時(shí)可以提取林地區(qū)域。其余地貌類型都滿足B2小于0.6，在此基礎(chǔ)上，當(dāng)B1大于0.35時(shí)可將冬小麥-花生和冬小麥-玉米大豆和其他作物和地貌類型做區(qū)分，根據(jù)B4大于0.20可以提取冬小麥-玉米大豆，B4小于0.20可以提取冬小麥-花生。其余作物和地貌類型滿足B1小于0.35，在此基礎(chǔ)上當(dāng)B3大于0.65時(shí)可將夏玉米和春玉米提取出來，此時(shí)當(dāng)B2大于0.45可提取春玉米，B2小于0.45可提取夏玉米。剩下兩種植被類型山地-夏玉米和雜草都滿足B3大于0.41且小于0.65，再根據(jù)B4大于0.20提取雜草，B4小于0.20提取山地-夏玉米。最后的地貌類型為非植被地表，在上述基礎(chǔ)上，B1大于0可提取居民地道路裸地，B1小于0可提取河流湖泊。

3 結(jié)果分析與精度驗(yàn)證

3.1 分類結(jié)果及分析

根據(jù)多時(shí)相遙感數(shù)據(jù)獲取的NDVI確定的最佳分類閾值，采用決策樹分類方法得到2019年和2020年9種地貌類型的分類結(jié)果如表2所示。

表2 2019年和2020年各地表類型分類占比 %

決策樹分類結(jié)果如圖5、圖6所示。

圖5 2019年決策樹分類結(jié)果

圖6 2020年決策樹分類結(jié)果

根據(jù)分類結(jié)果可知，與2019年相比，2020年林地面積增加了1.6%，造成略微增加的原因是2019年8月影像西北地區(qū)存在部分云，導(dǎo)致西北部分林地未分類(黑色區(qū)域)；河流湖泊下降0.2%，造成差異的原因是部分河流湖泊表面漂浮植被和藻類，使得少部分錯(cuò)分為植被；冬小麥-花生增加了4.2%，冬小麥-玉米大豆減少了3.8%，第一年雙熟制耕地大部分種植冬小麥-玉米大豆導(dǎo)致土壤中的磷元素過度消耗，為了避免連作對土壤造成傷害(土壤板結(jié)、土壤團(tuán)粒結(jié)構(gòu)破壞、土壤疏松度下降)，因此第二年雙熟制耕地大部分改種植冬小麥-花生；同理，春玉米增加了17.4%，夏玉米減少了19.6%也是為了降低對土壤的破壞；山地-夏玉米由于多種植在丘陵山區(qū)且密度較低，對土壤的傷害可忽略不計(jì)，所以種植面積基本無變化(增加了0.06%)；居民地道路裸地增加了2.8%，主要是城鎮(zhèn)化水平的提高導(dǎo)致居民地面積的增加；雜草減少了3.9%，雜草種類多密度低分散廣較難區(qū)分，每年的面積大致維持在10%左右。

3.2 精度驗(yàn)證

圖7 驗(yàn)證樣本點(diǎn)位置

分別在2019年8月和2020年8月(植被最茂盛期)進(jìn)行實(shí)地調(diào)查，再結(jié)合哨兵影像目視判讀等方式解譯出地表類型，以一個(gè)像元為基本單位建立驗(yàn)證樣本集，由于各地表類型所占面積差異過大，因此針對占比大的地貌類型選取較多驗(yàn)證樣本點(diǎn)，占比少的地貌類型選取較少驗(yàn)證樣本點(diǎn)(圖7)，以kappa系數(shù)、制圖精度和用戶精度對兩年的分類結(jié)果進(jìn)行驗(yàn)證。兩年的kappa系數(shù)分別為0.885 8和0.910 0，總體分類精度分別為90.76%和92.54%，9種植被類型的制圖精度和用戶精度如表3所示。根據(jù)kappa系數(shù)的判別原則：0～0.40(不含0.40)分類一致性差，0.40～0.80(不含0.80)分類中度一致，0.80～1分類高度一致，可知文中9種植被類型的遙感分類總體在精度上達(dá)到了高度一致。其中，各地貌類型樣本點(diǎn)被正確分類個(gè)數(shù)及占比如表4所示。

通過精度驗(yàn)證可知，林地、居民地道路裸地、河流湖泊面積大、地塊連片分布，有著明顯且易區(qū)分的NDVI值，面積估算精度很高，大致在95%左右；冬小麥-花生和冬小麥-玉米大豆和其他作物相比有著獨(dú)特的雙峰特性，面積估算精度較高，大致在90%～95%之間；春玉米、夏玉米、山地-夏玉米單熟制作物種植較分散，NDVI變化曲線一致，根據(jù)峰值特性進(jìn)行區(qū)分，面積估算精度一般，大致在80%左右；雜草由于其野生性，分布極其分散，很多雜草區(qū)域長寬小于10 m，小于哨兵2號最小像元面積，會與周圍地物形成混合像元，因此面積估算精度較差，大致在70%左右。

表3 2019年和2020年9種植被類型分類的制圖精度和用戶精度 %

表4 2019年和2020年各地貌類型驗(yàn)證樣本點(diǎn)個(gè)數(shù)及正確分類個(gè)數(shù)

4 結(jié)束語

文中研究了基于2019年所監(jiān)測的多時(shí)相NDVI曲線來構(gòu)建閾值，根據(jù)2019年和2020年4月、6月、8月、10月共8期影像，利用決策樹方法將耕地類、林草類以及非植被地表類共9種植被類型進(jìn)行了分類，最后根據(jù)驗(yàn)證樣本集利用混淆矩陣對分類結(jié)果進(jìn)行精度驗(yàn)證，結(jié)果表明：

1)該方法較好地提取了河北省靈壽縣2019年和2020年兩年9種不同的植被信息，其中包括3種不同品種的玉米，該技術(shù)方法不但可以對不同作物及植被類型進(jìn)行分類，而且可以較好地進(jìn)行精細(xì)分類(3種玉米的分類)，同時(shí)監(jiān)測作物的時(shí)空變化規(guī)律，可以對來年的作物種植趨勢進(jìn)行預(yù)測，在一定程度上為政府部門的行政決策提供依據(jù)。

2)采用哨兵2號多波段構(gòu)建NDVI值，并結(jié)合植被的多時(shí)相NDVI構(gòu)建決策樹進(jìn)行分類，與傳統(tǒng)的監(jiān)督分類和非監(jiān)督分類相比，該方法更加注重挖掘和分析不同植被的光譜信息差異，增加了植被類型可分性的信息量，分類識別的針對性更強(qiáng)，顯著提升了分類精度。

3)基于哨兵2號時(shí)序植被指數(shù)的作物精細(xì)分類及時(shí)空監(jiān)測上還有很多方面需要改進(jìn)。8月份的影像大多含有云層，云層的存在極大降低了分類精度；部分河流湖泊邊緣表面漂浮著植被，會將該部分水體錯(cuò)分為植被；玉米雜草分布不規(guī)則，混合像元較多，對面積提取的精度造成一定的影響。

針對上述問題今后可以將研究工作轉(zhuǎn)向哨兵1號影像的分類，利用哨兵1號雷達(dá)數(shù)據(jù)具有穿透性的特點(diǎn)，可以避免云層的影響，同時(shí)構(gòu)建雷達(dá)反射系數(shù)對植被類型進(jìn)行分類，將兩種分類結(jié)果進(jìn)行比較以提高整體分類精度。其次通過上述研究可知作物的種植成數(shù)(即面積比例)和地貌的形狀指數(shù)對分類結(jié)果影響較大，作物所種植的區(qū)域越集中越齊整，面積提取的精度也就越高，因此也可以將研究重心放在農(nóng)作物的地塊特征上。