基于中尺度光譜和時(shí)序物候特征提取南方丘陵山區(qū)茶園

馬超， 楊飛， 王學(xué)成

(1.中國(guó)科學(xué)院地理科學(xué)與資源研究所資源與環(huán)境信息系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100101； 2.中國(guó)科學(xué)院大學(xué)，北京 100049)

0 引言

中國(guó)是最早發(fā)現(xiàn)和利用茶樹的國(guó)家[1]。近些年來，隨著經(jīng)濟(jì)快速發(fā)展，茶園種植面積顯著增長(zhǎng)。茶園作為中國(guó)南方的特色經(jīng)濟(jì)林之一，在當(dāng)?shù)亟?jīng)濟(jì)發(fā)展中起著重要作用[2]。監(jiān)測(cè)南方丘陵山區(qū)茶園的分布和動(dòng)態(tài)變化，可為政府有關(guān)部門提供決策依據(jù)，對(duì)合理種植茶園和環(huán)境治理有重要意義。

南方丘陵山區(qū)茶園空間分布的統(tǒng)計(jì)方法有人工實(shí)地考察和遙感解譯2種。人工實(shí)地考察準(zhǔn)確度高，但是耗費(fèi)大量人力、物力和財(cái)力，不適合大面積大尺度的茶園種植區(qū)提取。遙感技術(shù)具有覆蓋范圍大、重訪周期短、應(yīng)用成本低等優(yōu)勢(shì)，能及時(shí)、準(zhǔn)確地掌握土地利用和土地覆被的現(xiàn)狀及動(dòng)態(tài)變化情況[3-8]，因此遙感解譯是實(shí)現(xiàn)大面積、大尺度茶園種植區(qū)提取的重要手段。然而，目前對(duì)于茶園提取的研究相對(duì)較少[9]，主要研究有： Rao等[10]利用歸一化植被指數(shù)(normalized difference vegetation index, NDVI)、簡(jiǎn)單比值(simple ratio, SR)植被指數(shù)和轉(zhuǎn)換型植被指數(shù)(transformed vegetation index, TVI)等構(gòu)建模型預(yù)測(cè)印度茶園產(chǎn)量，對(duì)比模型預(yù)測(cè)結(jié)果發(fā)現(xiàn)，基于TVI的產(chǎn)量模型與綠茶產(chǎn)量相關(guān)性最高； Dihkan等[11]采用支持向量機(jī)的方法提取安納圖里亞東北的茶園，發(fā)現(xiàn)運(yùn)用紋理分類結(jié)果的Kappa系數(shù)達(dá)到0.97； Dutta等[12]利用監(jiān)督和非監(jiān)督分類方法提取印度茶園，結(jié)果顯示，光譜相似性會(huì)造成茶樹與木本植被混淆； 徐偉燕等[13]基于資源三號(hào)衛(wèi)星數(shù)據(jù)，對(duì)比神經(jīng)網(wǎng)絡(luò)和決策樹分類方法的茶園提取結(jié)果，發(fā)現(xiàn)決策樹方法利用了光譜和紋理特征，可以提高茶園提取精度。

現(xiàn)有茶園提取研究多基于高空間分辨率影像，利用光譜特征和紋理特征進(jìn)行分類，影像數(shù)據(jù)成本高。而利用中、低空間分辨率影像進(jìn)行僅基于光譜特征的茶園提取時(shí)，由于茶園的光譜特征與灌叢、果園和林地有相似性，會(huì)產(chǎn)生錯(cuò)分現(xiàn)象。總體上看，目前茶園提取方法的指導(dǎo)性和實(shí)用性還有待提高，尚缺乏普適性和系統(tǒng)性的提取方法。針對(duì)以上問題，本文利用多時(shí)相，多源和中、低空間分辨率遙感影像數(shù)據(jù)，考慮光譜特征、植被指數(shù)和物候特征，提出一種系統(tǒng)、普適的方法，以期為相關(guān)工作者提供參考。

1 研究區(qū)概況及數(shù)據(jù)源

1.1 研究區(qū)概況

本文研究區(qū)為福建省漳州市和安溪縣。漳州市位于福建省東南部，下轄華安縣、平和縣、南靖縣、漳浦縣、龍海市和漳州市區(qū)等10個(gè)縣市，位置介于E117°～118°，N23.8°～25°之間，屬南亞熱帶季風(fēng)氣候區(qū)，年平均溫度為21 ℃，年降雨量1 500 mm。2016年漳州全市茶園面積46.9萬畝[注]1畝=666.67 m2。安溪縣的西南部與漳州市鄰接，隸屬泉州市，位于E117.6°～118.3°，N24.8°～25.4°之間，屬南、中亞熱帶海洋性季風(fēng)氣候區(qū)。安溪縣以農(nóng)業(yè)為主，主產(chǎn)茶葉，2014年安溪縣茶園總面積達(dá)60萬畝。研究區(qū)如圖1所示。影像為2015年Landsat8 OLI B5(R)，B4(G)，B3(B)波段組合影像。

圖1 研究區(qū)遙感影像Fig.1 Image of the study area

1.2 數(shù)據(jù)源

實(shí)驗(yàn)數(shù)據(jù)包括Landsat8 OLI數(shù)據(jù)(2015年的第290天)、2015年全年MODIS MOD13Q1數(shù)據(jù)、野外全球定位系統(tǒng)(global positioning system,GPS)實(shí)地采樣數(shù)據(jù)、Google Earth影像數(shù)據(jù)、土地利用數(shù)據(jù)、數(shù)字高程模型(digital elevation model，DEM)數(shù)據(jù)、研究區(qū)行政區(qū)劃數(shù)據(jù)以及漳州市和安溪縣茶園統(tǒng)計(jì)數(shù)據(jù)等，具體數(shù)據(jù)參數(shù)見表1。

表1 數(shù)據(jù)源介紹Tab.1 Data source

2 研究方法

2.1 茶園提取尺度與數(shù)據(jù)分辨率

根據(jù)漳州市和安溪縣茶園面積統(tǒng)計(jì)數(shù)據(jù)，茶園總種植面積超過700 km2，單塊茶園面積達(dá)幾十畝，甚至幾百畝。本文對(duì)Google Earth、Landsat8 OLI和MODIS 3種影像數(shù)據(jù)進(jìn)行對(duì)比(圖2)，發(fā)現(xiàn)研究區(qū)茶園面積能夠在中、高空間分辨率影像中占據(jù)多個(gè)像元。因此，使用Landsat8 OLI和MODIS數(shù)據(jù)等中、低空間分辨率影像可用于提取大面積茶園。

本文研究目的是提取茶園空間分布，雖然利用Landsat8 OLI的光譜特征能夠?qū)⒉鑸@與居民地、水體、農(nóng)田等區(qū)分開，但由于果園、部分林地與茶園的單一時(shí)相光譜特征相似性，茶園提取的精度較低。相比于Landsat8 OLI數(shù)據(jù)，盡管MODIS數(shù)據(jù)空間分辨率較低，但其獲取數(shù)據(jù)周期短，具有多時(shí)相特征，能夠獲取植被指數(shù)時(shí)間序列曲線，提取植被物候參數(shù)。對(duì)比茶園與果園、部分林地生長(zhǎng)季物候參數(shù)的差異對(duì)Landsat8 OLI數(shù)據(jù)分類結(jié)果進(jìn)一步提取，能夠改善單一時(shí)相光譜相似性的問題，從而提高分類精度。

(a) Google Earth影像數(shù)據(jù)(b) Landsat8 OLI 數(shù)據(jù)(c) MODIS數(shù)據(jù)

圖2不同空間分辨率數(shù)據(jù)對(duì)比

Fig.2Comparisonofdifferentspatialresolutiondata

2.2 技術(shù)路線

茶園種植區(qū)提取流程如圖3所示。

圖3 茶園種植區(qū)提取流程Fig.3 Flow chart of tea plantation extraction

首先基于MODIS增強(qiáng)植被指數(shù)(enhanced vegetation index，EVI)和NDVI時(shí)序曲線確定最合適的分類時(shí)間窗口，由時(shí)間窗口選擇合適的Landsat數(shù)據(jù)，對(duì)有云遮擋區(qū)域，選擇相近時(shí)期或前一年同時(shí)期無云影像替換； 通過野外GPS實(shí)地采樣數(shù)據(jù)與Google Earth影像選取采樣點(diǎn)，建立茶園解譯標(biāo)志庫，篩選光譜、地形等特征； 然后基于eCognition軟件平臺(tái)，使用面向?qū)ο蠓诸惙椒ê蜎Q策樹分類模型對(duì)研究區(qū)地物進(jìn)行分類，利用物候特征改善初步分類結(jié)果以提高茶園提取精度； 最后使用采樣數(shù)據(jù)對(duì)分類結(jié)果進(jìn)行精度評(píng)價(jià)。

2.3 時(shí)間窗口選擇

植被在不同生長(zhǎng)季節(jié)和生育期具有不同生理特征，多時(shí)相數(shù)據(jù)時(shí)間變化曲線能夠充分表現(xiàn)出同一植被類型在不同生長(zhǎng)期及不同植被類型在同一生長(zhǎng)期的差異[14-16]。利用MOD13Q1 NDVI/EVI數(shù)據(jù)產(chǎn)品提取2015年研究區(qū)7類土地利用類型(林地、果園、農(nóng)田、居民地、水體、茶園和其他)的MODIS植被指數(shù)時(shí)序曲線(圖4)。從圖4可以看出： ①居民地和水體的植被指數(shù)明顯低于植被； ②林地的植被指數(shù)最高； ③農(nóng)田的植被指數(shù)整體低于其他植被，有明顯峰值； ④總體上看，果園與茶園全年植被指數(shù)曲線相似，地物分類會(huì)造成果園與茶園的誤分； ⑤從茶園MODIS植被指數(shù)時(shí)序曲線上看，茶園的植被指數(shù)在7—8月出現(xiàn)峰值，在12月和1月出現(xiàn)最低值。由曲線可知10月份各地物植被指數(shù)區(qū)別較大，所以選擇Landsat OLI影像數(shù)據(jù)時(shí)間為2015年的第290天。

(a) MODIS-NDVI(b) MODIS-EVI

圖42015年MODIS-NDVI和MODIS-EVI時(shí)序曲線

Fig.4MODIS-NDVIandMODIS-EVItime-seriesprofilesin2015

2.4 面向?qū)ο蟮臎Q策樹分類

借助Google Earth影像和野外實(shí)地考察數(shù)據(jù)，隨機(jī)采集樣本點(diǎn)1 102個(gè)，建立研究區(qū)的遙感解譯標(biāo)志庫，包括土地利用類型7類： 林地、果園、農(nóng)田、居民地、水體、茶園和其他。利用面向?qū)ο蟮臎Q策樹方法進(jìn)行分類的步驟為： ①影像多尺度分割； ②光譜、地形特征選??； ③建立分類規(guī)則； ④決策樹分類。其中，影像分割尺度的選擇會(huì)直接影響最終分類結(jié)果[17]，通過多次試驗(yàn)，確定最優(yōu)分割尺度為30，形狀因子權(quán)重為0.1，緊致度因子權(quán)重為0.5。

2.5 物候特征選取與分析

植被物候監(jiān)測(cè)指通過各種方法和手段來監(jiān)測(cè)植物生長(zhǎng)節(jié)律的季節(jié)性和年際變化[18-20]。物候監(jiān)測(cè)主要利用植被指數(shù)時(shí)序曲線，常用的植被指數(shù)有NDVI和EVI[21-23]，考慮到MODIS-EVI數(shù)據(jù)產(chǎn)品已進(jìn)行了較好的大氣校正，且EVI計(jì)算公式中有對(duì)背景信息校正的參數(shù)，可避免基于比值的植被指數(shù)飽和問題的影響[24]，數(shù)據(jù)質(zhì)量要優(yōu)于其他植被指數(shù)。因此本文使用MODIS-EVI時(shí)序數(shù)據(jù)提取物候參數(shù)。

利用MODIS-EVI時(shí)序曲線提取的物候參數(shù)包含11個(gè)參數(shù)(圖5)，其中a和b分別表示生長(zhǎng)季的開始和結(jié)束點(diǎn)；c和d為振幅高度的80%；e為生長(zhǎng)季最大值；f為生長(zhǎng)季振幅；g表示生長(zhǎng)季長(zhǎng)度；h為生長(zhǎng)季期間EVI曲線與基準(zhǔn)值之間的積分值；ac斜率和bd斜率的絕對(duì)值分別表示生長(zhǎng)季的增長(zhǎng)率和下降率； 基準(zhǔn)值為曲線左側(cè)最小值和右側(cè)最小值的平均值[25-27]。

圖5 物候參數(shù)示意圖Fig.5 Phenological parameters

由圖5可知，物候參數(shù)a,b,c,d和h與均與振幅有相關(guān)性，因此本文選擇除此以外的6個(gè)物候參數(shù)用于茶園與農(nóng)田、林地等地類的區(qū)分(表2)。

表2 不同地物物候參數(shù)對(duì)比Tab.2 Comparison of phenological parameters

對(duì)比不同地類物候參數(shù)發(fā)現(xiàn)： ①茶園的生長(zhǎng)季最大值為5 285.677，遠(yuǎn)遠(yuǎn)小于林地、果園和其他地物； ②茶園生長(zhǎng)季振幅為2 314.807，與林地農(nóng)田區(qū)別較大； ③茶園生長(zhǎng)季增長(zhǎng)率在各地類中最小； ④茶園生長(zhǎng)季長(zhǎng)度及生長(zhǎng)季下降率與其他類別差異不明顯； ⑤茶園生長(zhǎng)季基準(zhǔn)值為2 970.887，與農(nóng)田區(qū)別較大。因此，本文選擇生長(zhǎng)季最大值、振幅、生長(zhǎng)季增長(zhǎng)率和基準(zhǔn)值4個(gè)物候參數(shù)作為茶園提取的物候指標(biāo)。

2.6 決策樹模型構(gòu)建與物候參數(shù)提取

在eCognition軟件平臺(tái)中使用決策樹方法對(duì)研究區(qū)Landsat8 OLI影像進(jìn)行分類。在決策樹模型中，使用紅光波段(R)、綠光波段(G)、藍(lán)光波段(B)、近紅外波段(NIR)、海藍(lán)波段(Coastal)、植被指數(shù)(NDVI/EVI)和亮度值(Brightness)構(gòu)造決策樹。具體參數(shù)范圍和相應(yīng)分類結(jié)果見圖6。

圖6決策樹分類方法構(gòu)建

Fig.6Constructionofdecisiontreeclassier

結(jié)果表明，利用光譜特征和植被指數(shù)等能夠?qū)崿F(xiàn)對(duì)居民地、水體、農(nóng)田和大部分林地分類。但是，由于果園和茶園的“異物同譜”的現(xiàn)象，導(dǎo)致大量果園錯(cuò)分為茶園。針對(duì)果園、部分林地和茶園的混合現(xiàn)象，對(duì)比分析茶園與果園、林地的物候參數(shù)，發(fā)現(xiàn)它們之間存在明顯差異。通過樣本點(diǎn)獲取茶園的物候參數(shù)范圍為： 98<生長(zhǎng)季增長(zhǎng)率<566，4 056<最大值<6 073，1 409<振幅<3 070，2 189<基準(zhǔn)值<3 509； 通過樣本點(diǎn)獲取茶園DEM范圍為263

3 結(jié)果與分析

3.1 茶園分布

本文利用面向?qū)ο蟮臎Q策樹分類模型得到茶園分類初步結(jié)果(圖7(a))，利用地物物候參數(shù)改善初步分類結(jié)果得到最終結(jié)果(圖7(b))。

(a) 未使用物候參數(shù)(b) 使用物候參數(shù)

圖7茶園分布提取結(jié)果

Fig.7Extractionresultofteaplantations

根據(jù)統(tǒng)計(jì)和調(diào)查數(shù)據(jù)，漳州市茶園主要分布在華安縣、平和縣、南靖縣和漳浦縣，安溪縣種植大面積茶樹。從圖7(a)可以看出，提取的茶園主要分布于安溪縣、漳州北部、漳州西部和東南部； 從圖7(b)可以看出，茶園主要分布于安溪縣、漳州北部、漳州西北部和東南部?？傮w上初步分類結(jié)果和最終分類結(jié)果與實(shí)際漳州市和安溪縣茶園分布一致。由圖7對(duì)比，茶園初步提取結(jié)果面積明顯多于最終分類結(jié)果，主要是茶園和果園、林地的錯(cuò)分，通過對(duì)比Google Earth影像與分類結(jié)果發(fā)現(xiàn)： 物候參數(shù)疊加決策樹分類結(jié)果能夠?qū)⒐麍@錯(cuò)分為茶園、林地錯(cuò)分為茶園的現(xiàn)象剔除； 同時(shí)沒有對(duì)茶園產(chǎn)生影響。

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

利用Google Earth影像隨機(jī)選取采樣驗(yàn)證點(diǎn)(共選取378個(gè))，通過構(gòu)造混淆矩陣得到漳州市和安溪縣地物分類精度評(píng)價(jià)(表3)。

表3 精度評(píng)價(jià)Tab.3 Accuracy assessment

研究區(qū)的茶園分類的總體精度為85.71%，Kappa系數(shù)為0.83，整體精度滿足分析需求。另外，從表中可以看出： ①果園和茶園分類混淆較多，表明兩者的光譜特征相似容易混淆。②茶園分類用戶精度為90.00%，生產(chǎn)者精度為83.72%，除了與果園和少量林地有錯(cuò)分現(xiàn)象外，整體分類結(jié)果滿足精度要求。

根據(jù)漳州市和安溪縣的統(tǒng)計(jì)數(shù)據(jù)，2016年漳州市茶園面積為312.67 km2，2015年安溪縣茶園面積達(dá)到400 km2。使用本文方法提取茶園面積為漳州市270.67 km2，安溪縣為453.33 km2。根據(jù)調(diào)查，漳州市的華安縣、平和縣、南靖縣和漳浦縣是茶園主要分布區(qū)，統(tǒng)計(jì)數(shù)據(jù)和分類結(jié)果見表4，與實(shí)際茶園種植面積相比有一定的誤差，但是總體精度超過80%，滿足空間分析與實(shí)際應(yīng)用需求。

表4 統(tǒng)計(jì)數(shù)據(jù)與分類結(jié)果對(duì)比Tab.4 Comparison of statistical data and classification data (km2)

4 結(jié)論

本文對(duì)南方丘陵山區(qū)茶園進(jìn)行提取，基于影像的光譜特征、植被指數(shù)和地形特征，構(gòu)造決策樹分類模型，并針對(duì)茶園與其他植被混淆情況，使用物候參數(shù)進(jìn)一步改善提取結(jié)果，研究發(fā)現(xiàn)：

1)使用MODIS-NDVI植被時(shí)序曲線可以觀測(cè)不同地物差異，發(fā)現(xiàn)漳州市及安溪縣的茶園提取時(shí)間窗口為第290天左右，茶園與其他地物有顯著區(qū)別，但由于“異物同譜”現(xiàn)象仍有果園、少部分林地與茶園之間的混淆。

2)決策樹分類模型注重對(duì)特征的分析和挖掘，使用光譜特征和植被指數(shù)能夠較好區(qū)分居民地、水體、林地和茶園，而且決策樹方法操作簡(jiǎn)單，精度較高。

3)使用MODIS-EVI時(shí)序數(shù)據(jù)提取物候參數(shù)，考慮了植被物候特征，改善了茶園分類結(jié)果，使得茶園總精度達(dá)85.71%，Kappa系數(shù)達(dá)0.83，分類結(jié)果滿足研究需要。

本文所提出的基于中尺度光譜和時(shí)序物候特征的茶園提取方法，適用于茶園提取及相關(guān)分析研究。未來可以進(jìn)一步考量不同分類模型和有關(guān)統(tǒng)計(jì)資料來實(shí)現(xiàn)更高精度提取，從而滿足土地利用和合理開發(fā)的更高要求。