柏 佳，孫 睿※，張赫林，王 巖，金志鳳

（1.北京師范大學(xué)地理科學(xué)學(xué)部，遙感科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京 100875；2.北京師范大學(xué)地理科學(xué)學(xué)部，北京市陸表遙感數(shù)據(jù)產(chǎn)品工程技術(shù)研究中心，北京 100875；3.中國科學(xué)院空天信息創(chuàng)新研究院/數(shù)字地球科學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100094；4.中國科學(xué)院大學(xué)，北京 100049；5.浙江省氣候中心，杭州 310002）

0 引 言

茶樹是中國重要的經(jīng)濟(jì)樹種，近年來茶園面積大幅增長[1-2]。茶樹的擴(kuò)張種植能夠促進(jìn)區(qū)域的經(jīng)濟(jì)發(fā)展，但是單一的樹種會(huì)降低生態(tài)系統(tǒng)的生物多樣性，加速水土流失、土壤酸化等生態(tài)環(huán)境問題[1,3-4]，及時(shí)掌握茶園的空間分布信息對茶園的生態(tài)管理、環(huán)境保護(hù)以及茶葉災(zāi)害的監(jiān)測預(yù)警具有重要意義。中國茶園大多分布在地形復(fù)雜的山區(qū)，地塊破碎，分布零散，形狀差異大、植被混雜且茶園所處環(huán)境長期受到云雨的影響，增加了茶園遙感提取的難度與不確定性，因此有必要研究遙感提取茶園的方法，來獲得較好的茶園識(shí)別精度。

近年來，遙感影像成為茶園提取的常用的數(shù)據(jù)來源[5]，例如，馬超等[6]利用30 m空間分辨率的Landsat8和250 m空間分辨率的MODIS數(shù)據(jù)對福建省安溪縣茶園進(jìn)行了提取，總體分類精度達(dá)到 85.71%；Zhu等[7]基于多時(shí)相 Sentinel-2的 10 m多光譜數(shù)據(jù)，基于隨機(jī)森林提取茶園，生產(chǎn)者精度和用戶精度分別達(dá)到96.57%和96.02%。但由于中國茶園大多分布在地形復(fù)雜的山區(qū)，地塊破碎，分布零散，形狀差異大，10～30 m空間分辨率的遙感數(shù)據(jù)不足以分辨小面積的茶樹種植區(qū)[8-10]。隨著遙感技術(shù)的發(fā)展，地物細(xì)節(jié)突出的更高空間分辨率遙感影像為茶園識(shí)別提供了新途徑[10-11]，其包含的紋理信息是解決茶園與其他植被光譜混淆的有效方法[12]。一些研究結(jié)合高空間分辨率遙感數(shù)據(jù)的紋理特征和光譜特征成功提取了茶園并取得較好的分類精度[13-15]。但是，在高分辨率遙感圖像上，農(nóng)作物與茶園的光譜特征和紋理特征十分接近，難以區(qū)分。歸一化差值植被指數(shù)（Normalized Difference Vegetation Index，NDVI）的時(shí)序信息可以反映茶園與農(nóng)作物的季節(jié)變化及物候期，通過比較茶園與農(nóng)作物 NDVI季節(jié)變化的相似性，可將部分農(nóng)作物與茶園區(qū)分開來。因此，有必要結(jié)合茶園與農(nóng)作物 NDVI的時(shí)序相似性提高二者的可分性[16-20]?；趧?dòng)態(tài)時(shí)間序列歸整算法（Dynamic Time Warping，DTW）計(jì)算的DTW距離能夠有效地度量兩個(gè)時(shí)序數(shù)據(jù)的相似性，目前已被證明是基于時(shí)序特征研究地表覆蓋分類的有效信息[21]。然而，高空間分辨率的衛(wèi)星重訪周期長，且茶樹種植區(qū)多云霧、多陰雨的特點(diǎn)造成時(shí)序數(shù)據(jù)的大量缺失[22-23]。目前，基于時(shí)空融合算法，將高空間分辨率數(shù)據(jù)的空間信息與中低空間分辨率的時(shí)間信息融合是解決時(shí)序數(shù)據(jù)難以獲取的有效手段[24]。

綜上所述，中低空間分辨率的遙感影像難以識(shí)別小面積的茶園，單一時(shí)相的高空間分辨率遙感影像的光譜信息和紋理信息無法很好地區(qū)分部分農(nóng)作物與茶園。表征NDVI時(shí)序相似性信息的DTW距離是提高農(nóng)作物與茶園可分性的重要特征，但受云霧的影響，完整的 NDVI時(shí)序數(shù)據(jù)難以獲取。針對以上問題，本文基于時(shí)空融合算法構(gòu)建時(shí)序NDVI，并結(jié)合高空間分辨率遙感影像的光譜特征和紋理特征研究茶園提取的方法，旨在及時(shí)掌握茶園分布，以期為茶園的管理及政策制定提供參考。

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

1.1 研究區(qū)概況

研究區(qū)王宅鎮(zhèn)隸屬浙江省金華市武義縣，位于武義縣中心，面積達(dá)101.7 km2，為武義縣第二大鎮(zhèn)。王宅鎮(zhèn)境內(nèi)四周丘陵起伏，西北地勢高峻，多峽谷丘陵，中部平坡盆地，土壤肥沃。氣候?yàn)閬啛釒Ъ撅L(fēng)氣候，氣候溫和，光照充足，雨量充沛，年平均溫度 16.7 ℃，年平均降水量1 450 mm。茶樹和農(nóng)作物為當(dāng)?shù)刂饕慕?jīng)濟(jì)作物，茶樹屬于常綠灌木，大多種植在丘陵地帶，農(nóng)作物類型主要包括水稻、玉米、馬鈴薯、油菜等，主要分布在平原地區(qū)。研究區(qū)如圖1所示。

1.2 遙感數(shù)據(jù)和預(yù)處理

所用數(shù)據(jù)包括 GF-1數(shù)據(jù)、Sentinel-2地表反射率數(shù)據(jù)和MODIS地表反射率產(chǎn)品（MCD43A4）。

為了降低其他植被對茶樹識(shí)別的影響，同時(shí)考慮到遙感數(shù)據(jù)的質(zhì)量，本文選取2017年冬季（2017年12月28日）覆蓋研究區(qū)的GF-1 2 m全色和8 m多光譜數(shù)據(jù)為主要數(shù)據(jù)源?；诘孛嬲{(diào)研得到的地面控制點(diǎn)（Ground Control Point，GCP）結(jié)合更高空間分辨率的影像（Google Earth）對GF-1影像進(jìn)行幾何精校正，輔以ASTER DEM數(shù)據(jù)對GF-1數(shù)據(jù)進(jìn)行正射校正，然后利用Gram-Schmidt融合方法將2m分辨率全色波段和8 m分辨率多光譜數(shù)據(jù)融合，得到2 m分辨率多光譜影像。

由于研究區(qū)農(nóng)作物類型眾多，物候復(fù)雜，僅基于幾期 NDVI數(shù)據(jù)來區(qū)分農(nóng)作物與茶園難以實(shí)現(xiàn)，因此需要更加完整的 NDVI時(shí)序數(shù)據(jù)來提高農(nóng)作物與茶園的可分性。然而，由于研究區(qū)全年多云多雨，完整的Sentinal-2時(shí)序NDVI難以獲取，因此，本文利用時(shí)間分辨率為1 d的MODIS反射率產(chǎn)品MCD43A4對Sentine-2 數(shù)據(jù)進(jìn)行時(shí)空融合得到完整的時(shí)序NDVI。

2017年研究區(qū)僅有9期Sentinel-2數(shù)據(jù)可用，可用數(shù)據(jù)日期分別為2017年2月28日、2017年4月29日、2017年7月18日、2017年9月16日、2017年10月1日、2017年10月31日、2017年12月10日、2017年12月20日和2017年12月2日。本文基于Google Earth Engine（GEE，https://earthengine.google.com/）平臺(tái)對Sentinel-2和MCD43A4數(shù)據(jù)進(jìn)行質(zhì)量控制后分別計(jì)算并提取了研究區(qū)2017年可用的9期Sentinel-2和2017年全年MODIS NDVI時(shí)序數(shù)據(jù)；利用 IDL平臺(tái)對 MODIS NDVI時(shí)序數(shù)據(jù)進(jìn)行插值和Savitzky-Golay（SG）濾波；最后基于貝葉斯最大熵（Bayesian Maximum Entropy，BME）時(shí)空統(tǒng)計(jì)的多源遙感數(shù)據(jù)融合方法[25-26]，融合MODIS和Sentinel-2的NDVI得到時(shí)空完整、時(shí)間分辨率5 d的10 m空間分辨率NDVI時(shí)序，為了避免空間統(tǒng)計(jì)估計(jì)中的空間平滑效應(yīng)，本研究只采用了時(shí)間協(xié)方差進(jìn)行融合，并對融合結(jié)果進(jìn)行SG濾波[27]。

2 研究方法

利用基于時(shí)空融合算法重建的Sentinel-2 NDVI時(shí)序數(shù)據(jù)，采用DTW算法計(jì)算茶園樣本與待分像元的DTW距離，用來表征茶園與待分像元的時(shí)相差異；結(jié)合GF-1的光譜特征和基于灰度共生矩陣（Gray Level Co-occurrence Matrix，GLCM）計(jì)算的紋理特征，選擇兩種數(shù)據(jù)組合方案（包括和不包括DTW距離），利用隨機(jī)森林算法提取茶園并進(jìn)行精度檢驗(yàn)，技術(shù)路線如圖2所示。

2.1 分類特征提取

2.1.1 光譜特征和紋理特征

光譜特征包括GF-1多光譜數(shù)據(jù)的光譜信息（包括近紅外、紅光、綠光和藍(lán)光波段）和基于光譜信息計(jì)算的單期NDVI。紋理特征由GF-1全色波段數(shù)據(jù)，以15為紋理提取窗口，基于GLCM計(jì)算包括均值（mean）、方差（variance）、同質(zhì)性（homogeneity）、對比度（contrast）、非相似性（dissimilarity）、熵（entropy）、角二階矩（angular second moment）和相關(guān)性（correlation）8個(gè)紋理特征[28-29]。

2.1.2 NDVI時(shí)序相似性

研究區(qū)平原地區(qū)存在大面積的冬種農(nóng)作物，其紋理和光譜特征與茶園非常相似，易造成二者的混分。為有效區(qū)分茶園與農(nóng)作物，提高茶園的分類精度，本文采用DTW算法[21,30]，計(jì)算待分像元與茶園樣本NDVI時(shí)序的相似性，該算法以兩個(gè)時(shí)列之間的距離來度量時(shí)序的相似性，兩個(gè)時(shí)序的相似性越高，則二者時(shí)序變化越相似，屬于同一地物類型的可能性越大。DTW算法的特點(diǎn)是可以自動(dòng)尋找目標(biāo)曲線和參考曲線的最佳對齊方式，該方式對應(yīng)兩條曲線距離最短，最短距離即為DTW距離，也就是兩個(gè)時(shí)序的相似度，距離越短，相似度越大[31-32]。

2.2 分類方法及精度評價(jià)

為了評估加入 NDVI時(shí)序特征對茶園提取的影響，選擇兩種特征組合方案進(jìn)行分類。方案 1：僅基于 GF-1的光譜特征和紋理特征；方案2：GF-1數(shù)據(jù)的光譜特征、紋理特征結(jié)合DTW距離特征。分類方法采用隨機(jī)森林分類方法。隨機(jī)森林是廣泛使用的基于集合的分類器，由若干個(gè)分類決策樹組成，每個(gè)決策樹為像素提供一個(gè)類別標(biāo)簽，最終獲得投票數(shù)最多的類別標(biāo)簽即為該像素的類別[33-34]。隨機(jī)森林在處理高維數(shù)據(jù)時(shí)具有靈活性和速度優(yōu)勢，并且可以解決過度擬合的問題[35-36]。

結(jié)合實(shí)地調(diào)研及更高分辨率的Google Earth圖像選取茶園（49 128個(gè)像元）與非茶園（552 311個(gè)像元）訓(xùn)練樣本進(jìn)行分類，樣本分布如圖1，分類工作基于Python平臺(tái)實(shí)現(xiàn)。

為了驗(yàn)證分類精度，結(jié)合Google Earth對研究區(qū)的茶園進(jìn)行目視判讀，并利用目視判讀結(jié)果隨機(jī)生成1 008個(gè)茶園和1 095個(gè)非茶園驗(yàn)證圖斑，驗(yàn)證圖斑分布如圖1所示?；诨煜仃?，采用準(zhǔn)確率、錯(cuò)誤率、精確率、召回率和F1分?jǐn)?shù)進(jìn)行茶園識(shí)別精度評價(jià)[37]。

3 結(jié)果與分析

3.1 Sentinel-2 NDVI的時(shí)空融合結(jié)果與NDVI時(shí)序相似性分析

研究區(qū)平原地區(qū)存在大面積的冬種農(nóng)作物，其紋理和光譜特征與茶園非常相似，容易造成二者的混分。為了有效區(qū)分茶園與農(nóng)作物，提高茶園的分類精度，本文基于時(shí)空融合獲得完整的NDVI時(shí)序數(shù)據(jù)并利用DTW算法計(jì)算茶園樣本與待分像元的DTW距離表征NDVI時(shí)序的相似性。圖3a顯示了2月8日NDVI的融合結(jié)果，該日研究區(qū)被云層完全覆蓋，無可用的 Sentinel-2數(shù)據(jù)。圖3b為2月28日實(shí)際Sentinel-2的NDVI圖像，時(shí)空融合后的NDVI圖像完整，較高值分布在北部和南部山區(qū)，較低值分布在中部城區(qū)，分布較為合理，地物細(xì)節(jié)突出，與2月28日實(shí)際Sentinel-2的NDVI相比，二者空間分布較為接近，說明基于MODIS和Sentinel-2 NDVI時(shí)序數(shù)據(jù)的時(shí)空融合結(jié)果可以用于茶園的提取。

基于DTW算法計(jì)算茶園樣本與待分像元NDVI時(shí)序的DTW距離如圖4所示。茶園樣本與茶園的DTW距離基本小于4，并主要分布在山區(qū)；由于常綠林的NDVI時(shí)序與茶園相似，研究區(qū)西北部的林地與茶園樣本的DTW距離基本小于4，南部地區(qū)的林地與茶園樣本的DTW距離在4～8之間；農(nóng)作物與茶園的NDVI時(shí)序差異較大，二者NDVI曲線的相似性較低，平原區(qū)農(nóng)作物與茶園樣本的DTW距離基本大于4；水體、建筑物和裸地等非植被的NDVI時(shí)序與茶園具有顯著差異，其DTW距離基本大于16。DTW距離的計(jì)算結(jié)果細(xì)節(jié)突出，有助于區(qū)分農(nóng)作物與茶園。

3.2 分類結(jié)果與精度驗(yàn)證

3.2.1 茶園提取的空間分布結(jié)果

基于兩種特征組合茶園的分類結(jié)果如圖5所示?？傮w而言，兩種特征組合提取的茶園空間分布較為一致。可以看出，研究區(qū)的大片茶園主要分布在西北和中部山區(qū)，中部平原地區(qū)分布較零散，面積較小。對于山區(qū)大片的茶園，二者的分類結(jié)果相似，均能較好地提取出茶園；在平原地區(qū)，添加DTW距離的茶園分類面積明顯小于未添加DTW距離的分類面積，這主要是由于將光譜和紋理特征與茶園相似的部分農(nóng)田與茶園區(qū)分開來，說明添加DTW距離能夠有效降低茶園與農(nóng)作物的混分，提高茶園的分類精度。

3.2.2 精度評價(jià)

基于混淆矩陣的精度評價(jià)如表1所示。添加DTW距離的準(zhǔn)確率為 96.91%，錯(cuò)誤率為3.09%，茶園提取的精確率為89.00%，召回率為83.09%，F(xiàn)1分?jǐn)?shù)為0.86；未添加DTW距離的分類結(jié)果的準(zhǔn)確率94.72%，錯(cuò)誤率為5.28%，茶園提取的精確率為73.09%，召回率為84.61%，F(xiàn)1分?jǐn)?shù)為0.78。添加DTW距離的準(zhǔn)確率較未添加DTW距離提高2.19個(gè)百分點(diǎn)，錯(cuò)誤率降低2.19個(gè)百分點(diǎn)，精確率提高15.91個(gè)百分點(diǎn)，召回率略微下降（1.52個(gè)百分點(diǎn)），但是兩種方案的召回率相差不大，F(xiàn)1分?jǐn)?shù)提高0.06，說明添加DTW距離的分類性能較好，茶園提取的精確率較高。

表1 不同特征組合的茶園分類結(jié)果精度評價(jià)Table 1 Accuracy assessment for extraction of tea plantations based on different feature combinations

加入DTW距離茶園的提取結(jié)果召回率略微下降，造成召回率略微降低的原因可能有：1）為了降低茶園種植造成的水土流失等問題，茶園內(nèi)部和周圍會(huì)種植一些較高樹木，GF-1的空間分辨率較高，導(dǎo)致其內(nèi)部和周圍存在這些樹木的陰影。陰影的NDVI時(shí)序的DTW距離與茶園存在差異，添加DTW距離后，會(huì)減少部分原本是茶園的識(shí)別；2）DTW 距離是基于融合后的 NDVI計(jì)算的，其空間分辨率是10 m，與GF-1的2 m分辨率相比，較粗空間分辨率的DTW距離的邊界鋸齒狀現(xiàn)象更加明顯，因此添加DTW距離特征時(shí)，茶園提取的部分邊界不及未添加DTW距離時(shí)完整和平滑，因此會(huì)造成茶園的召回率略微下降。

添加DTW距離的茶園分類的精確率顯著提高（提高15.91個(gè)百分點(diǎn)），說明未加入DTW距離的茶園分類面積明顯偏大，將非茶園的地類（如農(nóng)作物）錯(cuò)分為茶園，茶園的多分現(xiàn)象導(dǎo)致精確率比較低；而加入 DTW 距離后，由于茶園的NDVI時(shí)序與農(nóng)作物NDVI時(shí)序存在明顯差異，會(huì)顯著降低農(nóng)作物與茶園的混分，減小了茶園的多分現(xiàn)象，分類結(jié)果更加可靠。

4 討 論

本研究主要探討利用高空間分辨率遙感影像結(jié)合時(shí)序信息提取茶園的方法。馬超等[6]使用中空間分辨率的Landsat影像提取山區(qū)茶園，分類總體精度為85.71%。楊艷魁等[38]使用GF-2影像提取茶園，加入紋理信息后總體精度達(dá)到89.8%。徐偉燕等[13]使用資源三號(hào)衛(wèi)星數(shù)據(jù)，加入紋理信息作為分類特征，山區(qū)與平原地區(qū)的茶園提取精度分別達(dá)到92.97%和95%。本研究僅基于GF-1光譜與紋理信息茶園提取的準(zhǔn)確率達(dá)到94.72%。先前茶園提取的研究中，首先確定茶園與其他植被有差異的時(shí)間窗口，并以該時(shí)間窗口的影像作為時(shí)序特征。例如，馬超等[6]使用時(shí)序MODIS植被指數(shù)確定基于Landsat影像的茶園提取的最適窗口為10月份，選擇該月份Landsat影像進(jìn)行茶園分類，造成了茶園與林地、果園及其他植被類型的混分。趙曉晴等[39]僅使用三期Sentinel-2A數(shù)據(jù)作為時(shí)序特征，使用不同植被在 3期的差異構(gòu)建決策樹。本研究加入時(shí)序特征主要區(qū)分茶園與農(nóng)作物，而不同的農(nóng)作物生長季不同，僅使用某一個(gè)時(shí)間窗口的數(shù)據(jù)或較少的時(shí)序數(shù)據(jù)可能并不適用于區(qū)域內(nèi)所有的農(nóng)作物。因此，本研究基于Sentinel-2反射率數(shù)據(jù)與MCD43A4反射率產(chǎn)品進(jìn)行時(shí)空融合獲得完整的 NDVI時(shí)序數(shù)據(jù)，加入NDVI時(shí)序信息茶園提取的準(zhǔn)確率達(dá)到 96.91%，較未添加NDVI時(shí)序信息提高了2.19個(gè)百分點(diǎn)，由于茶園占比較小，因此總體精度提高不多。但是添加 NDVI時(shí)序信息茶園提取的精確率較未添加NDVI時(shí)序信息提高15.91個(gè)百分點(diǎn)，說明添加時(shí)序信息能夠有效提高茶園提取的精度。

但是本研究仍存在一些不確定性。首先，茶樹種植區(qū)多云霧、多陰雨，難以獲取高質(zhì)量的時(shí)序數(shù)據(jù)，為了彌補(bǔ)時(shí)序數(shù)據(jù)難以獲取的不足，本研究利用中低分辨率的MODIS時(shí)序數(shù)據(jù)對Sentinel-2的數(shù)據(jù)進(jìn)行了時(shí)序重建。在云量較大的時(shí)期，MODIS數(shù)據(jù)本身存在誤差且MODIS與Sentinel-2數(shù)據(jù)的空間分辨率差距較大，盡管MODIS高時(shí)間分辨率的特點(diǎn)基本能夠捕捉植被物候的差異，但是融合效果還是存在一些不確定性。在未來的研究中，需考慮影像質(zhì)量對融合精度的影響，嘗試?yán)每臻g分辨率較高的數(shù)據(jù)來融合時(shí)間分辨率較低的GF-1數(shù)據(jù)，進(jìn)一步提高茶園提取精度。其次，本研究僅使用GF-1數(shù)據(jù)為主要數(shù)據(jù)源，目前已有增加了能夠有效反映作物特有光譜特性的“紅邊”波段的GF-6數(shù)據(jù)。在未來的研究中，可以嘗試結(jié)合GF-1和GF-6來彌補(bǔ)GF-1時(shí)間分辨率較低的不足，及結(jié)合GF-6的“紅邊”波段，來提高茶園的提取精度。

5 結(jié) 論

本文通過對MODIS地表反射率產(chǎn)品和Sentinel-2反射率數(shù)據(jù)的時(shí)空融合，在獲得10 m空間分辨率的NDVI時(shí)序數(shù)據(jù)的基礎(chǔ)上，結(jié)合GF-1的光譜信息和紋理信息，利用隨機(jī)森林方法進(jìn)行了茶園分類，結(jié)果表明：

1）高時(shí)間、低空間分辨率的MODIS數(shù)據(jù)和高空間、時(shí)序不完整的Sentinel-2數(shù)據(jù)的融合，進(jìn)而得到完整的時(shí)序歸一化差值植被指數(shù)（Normalized Difference Vegetation Index，NDVI），及在此基礎(chǔ)上計(jì)算的動(dòng)態(tài)時(shí)間序列歸整距離（Dynamic Time Warping distance，DTW distance）能夠捕捉不同地物NDVI時(shí)序的差異，可用于茶園提取。

2）基于GF-1光譜、紋理信息結(jié)合DTW距離茶園分類結(jié)果的準(zhǔn)確率、錯(cuò)誤率、精確率、召回率和F1分?jǐn)?shù)分別為 96.91%、3.09%、89.00%、83.09%和 0.86，僅基于GF-1光譜和紋理信息茶園分類的準(zhǔn)確率、錯(cuò)誤率、精確率、召回率和F1分?jǐn)?shù)分別為94.72%、5.28%、73.09%、84.61%和 0.78。相對于僅基于光譜和紋理特征的分類結(jié)果，添加DTW距離的茶園分類精度有較大幅度提高，精確率和F1分?jǐn)?shù)分別提高了15.91個(gè)百分點(diǎn)和0.06，主要是由于時(shí)序信息的添加降低了茶園與平原地區(qū)農(nóng)作物的混分，表明高分辨率結(jié)合時(shí)序信息是提高茶園分類精度的有效手段。