祁向前,孫德浩,賈連星
(1.龍巖學(xué)院 資源工程學(xué)院,福建 龍巖 364012;2.中國礦業(yè)大學(xué) 環(huán)境與測繪學(xué)院,江蘇 徐州 221116;3. 黑龍江科技大學(xué) 礦業(yè)工程學(xué)院,哈爾濱 150022)
以往對(duì)農(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)行作物識(shí)別與監(jiān)測是一種行之有效的方法。例如Jakubauskas等[1]采用諧波算法對(duì)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é)合作物的物候特征對(duì)植被類型進(jìn)行分類。Conese等[7]對(duì)遙感數(shù)據(jù)進(jìn)行了單時(shí)相和多時(shí)相的研究,發(fā)現(xiàn)多時(shí)相數(shù)據(jù)能夠極大改善分類精度。賈樹海等[8]針對(duì)3個(gè)不同時(shí)期的遙感影像,結(jié)合不同作物的物候特征提取其NDVI植被指數(shù),采用監(jiān)督分類的方法,較好地提取了花生的種植面積并制圖。
哨兵2號(hào)是由兩顆極軌衛(wèi)星組成的星座,具有高分辨率多光譜成像的特點(diǎn),一顆衛(wèi)星的重訪周期為10 d,兩顆互補(bǔ),重訪周期為5 d,在紅邊范圍內(nèi)有較多波段,有利于對(duì)植被的監(jiān)測。國內(nèi)外均有學(xué)者利用單時(shí)相哨兵2號(hào)數(shù)據(jù)實(shí)現(xiàn)了作物識(shí)別[9-10]。NDVI作為一種植被指數(shù),被廣泛應(yīng)用于反映作物在不同季節(jié)、不同生育期表現(xiàn)出來的生理特征,也有效的應(yīng)用于土地利用覆蓋監(jiān)測、植被覆蓋密度評(píng)價(jià)、作物識(shí)別和作物產(chǎn)量預(yù)報(bào)等[11-13]。哨兵數(shù)據(jù)的波段最高分辨率為10 m,單純依靠目視解譯對(duì)作物進(jìn)行分類較為困難,因此構(gòu)建多時(shí)相NDVI植被指數(shù)來研究就能夠很好地提高分類精度[14-16]。
上述文獻(xiàn)主要是針對(duì)不同種類農(nóng)作物進(jìn)行分類或者是對(duì)單一農(nóng)作物進(jìn)行面積提取,涉及到的農(nóng)作物精細(xì)分類(例如夏玉米、春玉米、山地夏玉米等)及多種農(nóng)作物種植面積變化研究還很少。文中以河北省石家莊市靈壽縣的多種作物及地貌類型為例,采用多時(shí)相決策樹分類法,利用所建立的作物及地貌樣本點(diǎn)不同時(shí)期內(nèi)監(jiān)測到的NDVI值建立時(shí)序NDVI曲線,并在具有明顯區(qū)分特征的月份確定閾值,對(duì)研究區(qū)域內(nèi)的多種作物及地貌類型進(jìn)行分類研究。
靈壽縣(113°45′E~114°28′E,38°16′N~38°48′N)位于河北省西部,太行山東側(cè),在省會(huì)石家莊西北35 km處,東南與正定縣毗連,南與鹿泉區(qū)隔滹沱河相望,西同平山縣接壤,西北隅與平山縣、阜平縣及山西省五臺(tái)縣接境(圖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ū)域
根據(jù)作物的物候特征了解到作物有明顯植被指數(shù)區(qū)分度的月份為4月、6月、8月和10月,后續(xù)用所監(jiān)測樣本點(diǎn)的時(shí)序植被指數(shù)變化分析也驗(yàn)證了該想法,并且考慮到云量會(huì)對(duì)影像分類精度造成影響,因此通過綜合比較各期影像云量覆蓋率,分別選取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號(hào)由哨兵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ā)送下載請(qǐng)求,過半個(gè)至一個(gè)工作日即可下載。
表1 影像清單及各農(nóng)作物發(fā)育期
從歐空局下載下來的數(shù)據(jù)是L1C級(jí),無需進(jìn)行輻射定標(biāo),經(jīng)過大氣校正后進(jìn)行重采樣(10 m),進(jìn)行波段融合,剔除掉不需要的海岸/氣溶膠波段(band.1)和水蒸氣波段(band.9)。經(jīng)過裁剪和空間配準(zhǔn),得到處理后的影像數(shù)據(jù)。
圖2 技術(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)證樣本集,對(duì)當(dāng)年的作物及其他地貌分類通過混淆矩陣進(jìn)行精度評(píng)價(jià)來驗(yàn)證分類的準(zhǔn)確性。
(1)
式中:BNIR為哨兵2號(hào)近紅外波段(第八波段)的地表反射率;BRED為哨兵2號(hào)紅光波段(第四波段)的地表反射率。
圖3 樣本點(diǎn)位置
本研究區(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ù)并不會(huì)出現(xiàn)大降幅,其中林地常年較高,雜草多出現(xiàn)在道路兩旁及山區(qū),植被指數(shù)在8月到達(dá)峰值時(shí)低于作物的峰值(圖4)。
圖4 農(nóng)作物及林草類NDVI時(shí)序曲線
在構(gòu)建的NDVI時(shí)序曲線中,非植被地表由于全年低植被指數(shù)因此有著極高的識(shí)別度;雙熟制耕地作物及林地特征差異明顯,識(shí)別度較高;3種玉米具有相同的植被指數(shù)變化趨勢,在起伏及峰值上有差異,也較容易識(shí)別;雜草種類多樣,密度低,分布不集中,不能獲得全部種類雜草的植被指數(shù)變化情況,只能得到一個(gè)大致變化范圍,識(shí)別度一般。綜上所述,利用NDVI植被指數(shù)來對(duì)作物及其他地表類型進(jìn)行識(shí)別分類具有較好的可行性。
決策樹作為一種分類方法,具有靈活、直觀、運(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)建決策樹對(duì)作物及地貌進(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可提取河流湖泊。
根據(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)致土壤中的磷元素過度消耗,為了避免連作對(duì)土壤造成傷害(土壤板結(jié)、土壤團(tuán)粒結(jié)構(gòu)破壞、土壤疏松度下降),因此第二年雙熟制耕地大部分改種植冬小麥-花生;同理,春玉米增加了17.4%,夏玉米減少了19.6%也是為了降低對(duì)土壤的破壞;山地-夏玉米由于多種植在丘陵山區(qū)且密度較低,對(duì)土壤的傷害可忽略不計(jì),所以種植面積基本無變化(增加了0.06%);居民地道路裸地增加了2.8%,主要是城鎮(zhèn)化水平的提高導(dǎo)致居民地面積的增加;雜草減少了3.9%,雜草種類多密度低分散廣較難區(qū)分,每年的面積大致維持在10%左右。
圖7 驗(yàn)證樣本點(diǎn)位置
分別在2019年8月和2020年8月(植被最茂盛期)進(jìn)行實(shí)地調(diào)查,再結(jié)合哨兵影像目視判讀等方式解譯出地表類型,以一個(gè)像元為基本單位建立驗(yàn)證樣本集,由于各地表類型所占面積差異過大,因此針對(duì)占比大的地貌類型選取較多驗(yàn)證樣本點(diǎn),占比少的地貌類型選取較少驗(yàn)證樣本點(diǎn)(圖7),以kappa系數(shù)、制圖精度和用戶精度對(duì)兩年的分類結(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號(hào)最小像元面積,會(huì)與周圍地物形成混合像元,因此面積估算精度較差,大致在70%左右。
表3 2019年和2020年9種植被類型分類的制圖精度和用戶精度 %
表4 2019年和2020年各地貌類型驗(yàn)證樣本點(diǎn)個(gè)數(shù)及正確分類個(gè)數(shù)
文中研究了基于2019年所監(jiān)測的多時(shí)相NDVI曲線來構(gòu)建閾值,根據(jù)2019年和2020年4月、6月、8月、10月共8期影像,利用決策樹方法將耕地類、林草類以及非植被地表類共9種植被類型進(jìn)行了分類,最后根據(jù)驗(yàn)證樣本集利用混淆矩陣對(duì)分類結(jié)果進(jìn)行精度驗(yàn)證,結(jié)果表明:
1)該方法較好地提取了河北省靈壽縣2019年和2020年兩年9種不同的植被信息,其中包括3種不同品種的玉米,該技術(shù)方法不但可以對(duì)不同作物及植被類型進(jìn)行分類,而且可以較好地進(jìn)行精細(xì)分類(3種玉米的分類),同時(shí)監(jiān)測作物的時(shí)空變化規(guī)律,可以對(duì)來年的作物種植趨勢進(jìn)行預(yù)測,在一定程度上為政府部門的行政決策提供依據(jù)。
2)采用哨兵2號(hào)多波段構(gòu)建NDVI值,并結(jié)合植被的多時(shí)相NDVI構(gòu)建決策樹進(jìn)行分類,與傳統(tǒng)的監(jiān)督分類和非監(jiān)督分類相比,該方法更加注重挖掘和分析不同植被的光譜信息差異,增加了植被類型可分性的信息量,分類識(shí)別的針對(duì)性更強(qiáng),顯著提升了分類精度。
3)基于哨兵2號(hào)時(shí)序植被指數(shù)的作物精細(xì)分類及時(shí)空監(jiān)測上還有很多方面需要改進(jìn)。8月份的影像大多含有云層,云層的存在極大降低了分類精度;部分河流湖泊邊緣表面漂浮著植被,會(huì)將該部分水體錯(cuò)分為植被;玉米雜草分布不規(guī)則,混合像元較多,對(duì)面積提取的精度造成一定的影響。
針對(duì)上述問題今后可以將研究工作轉(zhuǎn)向哨兵1號(hào)影像的分類,利用哨兵1號(hào)雷達(dá)數(shù)據(jù)具有穿透性的特點(diǎn),可以避免云層的影響,同時(shí)構(gòu)建雷達(dá)反射系數(shù)對(duì)植被類型進(jìn)行分類,將兩種分類結(jié)果進(jìn)行比較以提高整體分類精度。其次通過上述研究可知作物的種植成數(shù)(即面積比例)和地貌的形狀指數(shù)對(duì)分類結(jié)果影響較大,作物所種植的區(qū)域越集中越齊整,面積提取的精度也就越高,因此也可以將研究重心放在農(nóng)作物的地塊特征上。