基于GF-1和Sentinel-2時序數(shù)據(jù)的茶園識別

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

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

0 引 言

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

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

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

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

1.1 研究區(qū)概況

研究區(qū)王宅鎮(zhèn)隸屬浙江省金華市武義縣，位于武義縣中心，面積達101.7 km2，為武義縣第二大鎮(zhèn)。王宅鎮(zhèn)境內(nèi)四周丘陵起伏，西北地勢高峻，多峽谷丘陵，中部平坡盆地，土壤肥沃。氣候為亞熱帶季風(fēng)氣候，氣候溫和，光照充足，雨量充沛，年平均溫度 16.7 ℃，年平均降水量1 450 mm。茶樹和農(nóng)作物為當(dāng)?shù)刂饕慕?jīng)濟作物，茶樹屬于常綠灌木，大多種植在丘陵地帶，農(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ù)據(jù)的質(zhì)量，本文選取2017年冬季（2017年12月28日）覆蓋研究區(qū)的GF-1 2 m全色和8 m多光譜數(shù)據(jù)為主要數(shù)據(jù)源?；诘孛嬲{(diào)研得到的地面控制點（Ground Control Point，GCP）結(jié)合更高空間分辨率的影像（Google Earth）對GF-1影像進行幾何精校正，輔以ASTER DEM數(shù)據(jù)對GF-1數(shù)據(jù)進行正射校正，然后利用Gram-Schmidt融合方法將2m分辨率全色波段和8 m分辨率多光譜數(shù)據(jù)融合，得到2 m分辨率多光譜影像。

由于研究區(qū)農(nóng)作物類型眾多，物候復(fù)雜，僅基于幾期 NDVI數(shù)據(jù)來區(qū)分農(nóng)作物與茶園難以實現(xiàn)，因此需要更加完整的 NDVI時序數(shù)據(jù)來提高農(nóng)作物與茶園的可分性。然而，由于研究區(qū)全年多云多雨，完整的Sentinal-2時序NDVI難以獲取，因此，本文利用時間分辨率為1 d的MODIS反射率產(chǎn)品MCD43A4對Sentine-2 數(shù)據(jù)進行時空融合得到完整的時序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/）平臺對Sentinel-2和MCD43A4數(shù)據(jù)進行質(zhì)量控制后分別計算并提取了研究區(qū)2017年可用的9期Sentinel-2和2017年全年MODIS NDVI時序數(shù)據(jù)；利用 IDL平臺對 MODIS NDVI時序數(shù)據(jù)進行插值和Savitzky-Golay（SG）濾波；最后基于貝葉斯最大熵（Bayesian Maximum Entropy，BME）時空統(tǒng)計的多源遙感數(shù)據(jù)融合方法[25-26]，融合MODIS和Sentinel-2的NDVI得到時空完整、時間分辨率5 d的10 m空間分辨率NDVI時序，為了避免空間統(tǒng)計估計中的空間平滑效應(yīng)，本研究只采用了時間協(xié)方差進行融合，并對融合結(jié)果進行SG濾波[27]。

2 研究方法

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

2.1 分類特征提取

2.1.1 光譜特征和紋理特征

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

2.1.2 NDVI時序相似性

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

2.2 分類方法及精度評價

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

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

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

3 結(jié)果與分析

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

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

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

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

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

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

3.2.2 精度評價

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

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

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

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

4 討 論

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

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

5 結(jié) 論

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

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

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