廖廓, 聶磊, 楊澤宇, 張紅艷, 王艷杰,彭繼達, 黨皓飛, 冷偉
(1.福建省氣象科學研究所,福州 350008; 2.武夷山國家氣候觀象臺,武夷山 354300; 3.武漢珈和科技有限公司,武漢 430200; 4.浙江萬維空間信息技術有限公司,南京 210012)
福建省武夷山市是中國著名茶鄉,也是世界紅茶和烏龍茶的發源地,“萬里茶道”的起點。茶葉是當地重要的產業支柱,受到政府單位的高度重視。如何對茶園進行規范管理、及時防治病蟲害、提高茶葉產量,是亟待解決的問題。及時掌握茶樹種植分布是監測茶園情況的基礎。近年來,衛星遙感技術在自然資源管理,農作物估產、長勢監測等工作中已成為主要手段,較傳統的人工野外測繪方式具有明顯的質量和成本優勢,得到了日益廣泛的應用。
目前,很多學者借助遙感影像豐富的光譜或空間信息對茶園分布進行提取。已有采用遙感影像進行茶園種植分布提取的研究主要分為基于像元的分類與基于對象的分類。基于像元的分類方法以像元為分類基本單位,根據像元的光譜、紋理、空間等特征,采用監督分類或決策樹分類方法判斷該像元是否為茶園區域。這種傳統的遙感影像分析方法主要被應用于中等分辨率遙感影像(如GF-1 WFV,Landsat TM,Landsat OLI,HJ-1A/B,Sentinel-2)的茶園分布提取,分類算法以支持向量機、隨機森林等淺層機器學習方法為主[1-9]。這類基于像元的方法很少考慮到相鄰像元間的空間結構信息,未能充分利用影像空間紋理信息,分類結果經常出現“椒鹽”噪聲[10]。在高分辨率遙感影像上,地物的空間紋理信息更加豐富,其細節特征也更加突出,使得類內方差增大、類間方差減小,基于像素的分類方法往往會引起大量的地物錯分。
常用于高分辨率影像的提取方法是基于對象的分類方法,首先對遙感影像進行分割,提取影像特征,構建對象,再以對象為基本分類單位,判斷對象是否為茶園[11-15]。相比于基于像元的分類技術,一定程度上能夠克服“椒鹽”現象等缺陷[11]。但該類方法不能充分學習同類對象特征的結構信息和規律,因而最后得到的分類精度也不樂觀,同時該方法依賴于分割方法以及分割參數的選取,普適性很差[16]。茶樹具有特殊的種植、管理方式,因此茶園在高分辨率影像上有特殊的特征,在高分辨率影像上茶園應是由多種地物構成的語義場景,而非一個簡單的土地覆蓋類別[17]。近年來,越來越多的學者嘗試從場景尺度對茶園進行提取[18-19],并取得了不錯的效果。在這些研究中,大多數基于對象的影像分類方法,都采用人工特征建模,即通常采用諸如顏色、紋理、幾何等低層影像特征,這些低層特征不能真正表示影像中感興趣的類別。在基于場景的研究中,最優視覺單詞數和主題個數的確定,也需要人工過多的干涉。
近年來,隨著深度學習技術的發展,越來越多的卷積神經網絡(convolutional neural network,CNN)方法被成功應用高分辨率遙感影像的茶園分類研究。陳厚坤等[20]使用0.26 m Google影像,基于預訓練的VGG-16網絡結構,采用遷移學習方法訓練茶園樣本進行茶園分類研究; Tang等[21]提出一種基于對象的CNN,采用0.5 m分辨率的Google影像訓練模型,使用遷移學習的方式,將模型應用到1 m GF-2影像的茶園分類; Jamil等[22]搭建6層CNN,從高分辨率航空影像中提取概括性高且魯棒性強的影像特征,采用隨機森林對特征進行分類以提取茶園。
整體上,目前茶園提取方法的普適性以及實用性還有待提高,尚缺乏普適性和系統性的高精度提取方法。基于多源遙感和時間序列數據,集成智能算法對茶園進行大尺度遙感監測將成為必然的發展趨勢[23]。本文針對單一影像源茶園難提取問題,綜合多光譜影像的光譜信息以及超高分辨率影像的紋理特征,提出一種基于多源高分辨率衛星影像和多維卷積神經網絡(multidimensional multi-source convolutional neural networks, MM-CNN)的茶園分類方法。以一維和二維CNN(1D-CNN,2D-CNN)為基礎,采用10 m分辨率Sentinel-2影像以及0.5 m分辨率Google影像,通過建立2種模型,分別提取茶園及疑似區域,融合2個模型結果,得到最終茶園分布。
本文研究區域為福建省武夷山市新田鎮,經緯度范圍為E117.84°~118.15°, N27.45°~27.64°。地形地貌復雜(圖1)。區域氣候屬中亞熱帶季風濕潤氣候,四季分明,雨量豐沛。該區域為茶園種植大區,區域內地物類型豐富,包括林地、茶園、耕地、水域、草地、居民地、道路,對研究茶園與林地、耕地等背景地物的分類較適宜。為確保茶產業可持續發展,近年來,武夷山市持續對茶山進行綜合整治,新增茶園較少[24]。

圖1 研究區域Google影像、野外調查圖斑及地面真實標簽制作區域
本研究使用的數據有遙感數據以及野外實地調查數據。其中遙感數據包括Google影像、Sentinel-2影像,空間分辨率分別為0.5 m和10 m。受當地云雨天氣影響,從歐洲航天局哥白尼數據中心(https: //scihub.copernicus.eu/)僅查詢、獲取到2019年全年5期無云Sentinel-2影像,分別是20190127,20190919,20191019,20191113和20191213,其中20190127期影像有薄云。區域內Google影像由2018—2020年多期影像合成,不同時期影像色調不一致導致研究區Google影像具有明顯的鑲嵌痕跡(圖1)。為將本文方法推廣至中高分辨率衛星影像(如Sentinel-2,Landsat8,GF-1,GF-6等)的茶園分布提取,考慮不同衛星影像波段差異,本文僅選用Sentinel-2影像藍光(459.4~525.4 nm)、綠光(541.8~577.8 nm)、紅光(649.1~680.1 nm)和近紅外(779.8~885.8 nm)4個中高分辨率衛星影像都具有的波段進行研究。原始L1C級Sentinel-2影像經Sen2Cor軟件輻射定標、大氣校正后再進行波段組合、裁剪、拼接得到研究區Sentinel-2影像。
研究所使用的野外調查數據是于2019年5月在該研究區域內通過手持GPS儀器采集地物樣本及控制點信息獲取的。地物采集內容主要包含茶園、耕地、林地、草地、水域、園林(非茶園園林)、設施用地、其他共8種地物類型,通過內業編輯生成77個圖斑,調查總面積357 454 m2,其中茶園圖斑25個(占比32%),面積95 293 m2(占比27%)。根據實調結果,在實調區域內選取4個矩形區域制作茶園地面真實標簽,如圖1所示。
茶樹屬于灌木或小喬木,四季常綠,與灌叢、喬木等植被光譜特征相似。受人為經營管理影響,茶樹新梢一年中在3—10月生長多次,11月—次年2月則處于休眠期。可多次采摘,在3—4月采摘春茶,9月采摘秋茶。茶樹修剪在采摘后進行,一般集中在5—6月。
結合不同地物物候特征以及遙感影像可獲取性,基于野外實調數據,采用2019年1月、9月、10月、11月以及12月中旬Sentinel-2遙感影像,分別進行不同地物光譜分析。通過提取不同地物在不同波段的光譜反射率,繪制不同時期不同地物類型的均值光譜曲線圖(圖2)以及不同波段光譜反射率隨時間變化曲線(圖3),對比分析茶園、其他地物光譜特征及其差異。冬季1月,茶樹處于休眠期,在可見光波段,茶園光譜反射率明顯低于除林地外的其他地物; 在近紅外波段茶園光譜反射率僅低于設施用地,高于其他地物,該時期茶園與其他地物光譜差異較明顯。夏末初秋9月,各植被類型都表現出生長旺盛的植被特征,沒有清晰的光譜差異,但部分秋茶采摘會導致茶園植被特征減弱,從而改變其光譜特征。初秋10月,由于各植被生長物候的差異,該時期茶樹生長減緩,茶園與其他植被有一定的差異,在可見光波段光譜反射率低于除林地外的其他植被,在近紅外波段高于其他植被,但與草地以及耕地較接近。深秋11月,由于作物的收割,該時期茶園在近紅外波段光譜反射率與耕地有了較明顯差異,但與草地、林地差異很小。初冬12月,茶樹生長處于休眠期,不同植被物候差異進一步增強,該時期茶園僅在近紅外波段上與林地有較小的差異。茶園、林地、園林3種植被一年中在不同波段的反射率變化趨勢較一致(圖3),但在藍光、綠光波段的反射率差異較小,在紅光以及近紅外波段反射率差異較大,其中茶園與園林紅光波段反射率在12月差異最大,茶園與林地近紅外波段光譜反射率在10月差異最大。

(a) 1月(b) 9月(c) 10月(d) 11月(e) 12月

圖2 實調區地物在不同月份光譜曲線圖

(a) 藍光波段(b) 綠光波段(c) 紅光波段(d) 近紅外波段

圖3 實調區地物光譜反射率時間變化曲線
總體上,研究區內不同植被在不同時期均表現出較明顯的季節光譜特征,茶園與其他植被光譜特征在9月影像上差異很微弱,在1月差異最明顯。因此1月是茶園分類提取的最佳時期,但該時期影像質量較差。綜合考慮影像質量與各地物類型光譜特征差異,最終選取10月以及12月影像提取茶園分布,該時期茶樹處于生長休止期,且茶園光譜不受采摘、修剪等農事活動干擾。
茶樹常見種植方式為行播,茶園在0.5 m空間分辨率的Google影像上呈現出條狀紋理特征,在10 m空間分辨率的Sentinel-2影像上無條紋狀紋理特征。由于研究區Google影像為不同時期影像拼接而成,不同時期影像成像條件差異以及地形起伏會帶來的不同影像畸變,以及新老茶園特征不一導致茶園在影像上紋理特征多樣化。研究區Google影像茶園典型紋理特征如圖4所示。圖4(a)是典型老茶園條狀紋理特征,該特征與部分規模化種植的經濟作物(如花卉、竹蓀等)相似; 圖4(b)植被特征較弱,但具有明顯的條狀紋理特征,為典型新茶園紋理特征,該種特征與休耕狀態的梯田相似; 圖4(c)所示的茶園部分呈現條狀紋理特征,部分不具有條狀紋理,這是不同時期影像拼接邊緣羽化導致的,茶園紋理平滑后的特征與部分耕地、水域紋理特征相似。

(a) 老茶園紋理(b) 新茶園紋理(c) 模糊紋理
1D-CNN常用于時間序列數據的處理[25]。本文構建的1D-CNN模型網絡結構如圖5所示,由4個5×5卷積層(“same”填充,線性激活函數(ReLU))、4個2×2最大池化層以及全連接層組成。池化層連接在卷積層后,步長設置為2,起降低向量大小、減小權值參數作用,防止計算過程中出現過擬合現象。進入第一個全連接層前,對輸入進行扁平化處理,第一個全連接層設置512個神經單元,采用ReLU線性激活函數; 最后一個全連接層神經元個數設置為2,采用softmax激活函數輸出最終分類結果。

圖5 1D-CNN網絡結構
U-Net網絡[26]是一種2D-CNN網絡結構,本文使用該網絡作為高分辨率遙感影像茶園分割模型。如圖6所示,網絡由一條收縮路徑(左側)與一條擴張路徑(右側)組成。收縮路徑是一個典型的卷積網絡結構,由2個3×3卷積層(無填充卷積,ReLU函數)以及一個2×2最大池化重復組成。擴張路徑中的每個步驟首先對特征圖進行上采樣,然后進行2×2卷積(向上卷積),以將特征圖通道的數量減半,并與從相應收縮路徑中裁剪特征圖進行級聯,再進行2個3×3卷積,每個卷積后是ReLU激活函數。由于每次卷積后邊緣像素都會丟失,因此對收縮路徑中的特征圖進行裁剪是必要的。在最后一層,使用1×1卷積將每個64分量特征向量映射到所需的類別數。網絡共23個卷積層。

圖6 U-Net網絡結構
基于MM-CNN的茶園分類方法流程如圖7所示。從多期Sentinel-2茶園樣本中獲取茶樹物候信息,確定不同時相的解譯特征,人工繪制樣本,應用于1D-CNN模型中,完成Sentinel-2茶園提取模型的構建; 將Google影像茶園樣本紋理信息應用于U-Net分割模型,構建Google茶園提取模型; 多期Sentinel-2影像、Google影像分別輸入對應茶園提取模型,預測Sentinel-2影像、Google影像茶園分布。對Sentinel-2影像、Google影像茶園分類結果進行融合、后處理得到最終茶園分布結果。

圖7 本研究茶園分布提取流程
不同影像分類結果融合過程為: 首先對Sentinel-2影像、Google影像茶園分類結果進行后處理; 其次使用最近鄰重采樣將Sentinel-2影像分類后處理結果(10 m)重采樣至0.5 m; 然后對Sentinel-2影像分類后處理結果(0.5 m)與Google影像分類后處理結果(0.5 m)求交,得到融合后結果; 最后再次對融合結果進行后處理,得到最終分類結果。其中,分類結果采用小圖斑去除、空洞填充、邊緣平滑的等后處理方式優化(圖8),其中Sentinel-2影像分類結果最小碎斑點大小為5、最小空洞為2; Google影像分類結果最小碎斑點大小為200、最小空洞大小為100; 融合結果最小碎斑點大小為2 000、最小空洞為50,采用PAEK指數多項式平滑算法[27]對融合結果邊緣進行平滑。

(a) 處理前(b) 處理后
本文采用Python編程語言基于TensorFlow、Keras深度學習庫構建1D-CNN模型以及U-Net框架。參照楊澤宇等[28]的研究,1D-CNN模型學習率設置為1e-4,分塊大小為1 000,訓練樣本進行標準差歸一化預處理,模型循環次數設定為1 000。U-Net模型優化器采用自適應矩估計算法(Adam)來動態調整迭代參數的學習率。模型初始學習率均設置為0.001,選用類平衡交叉熵算法作為損失函數來對模型參數進行優化。訓練樣本在輸入模型訓練之前,對樣本作去中心化處理,即所有樣本減去樣本均值。U-Net模型循環次數設置為30。
根據武夷山市茶園實調數據建立不同影像解譯標志,繪制茶園與其他地物標準樣本。基于2019年10月19日、2019年12月13日Sentinel-2影像,分別繪制14 521個像素(茶園樣本4 523個像素)、16 542個像素(茶園樣本4 987個像素)的1D-CNN模型訓練樣本。基于Google影像繪制20個9 000像素×9 000像素大小影像標簽(包括茶園與非茶園),通過隨機裁剪、樣本擴增得到18 000對(影像+標簽)512×512像素大小的UNet模型訓練樣本。考慮到Google影像的色調差異(圖1),并最大化保留影像紋理特征,樣本擴增方式為翻轉、旋轉、增強對比度的隨機組合。所有模型樣本均按8∶2的比例隨機抽樣生成訓練集與驗證集樣本。
為檢驗不同模型的預測效果,我們將預測結果與實地調研結果進行對比,基于混淆矩陣計算F1值、總體精度(overall accuracy, OA)以及Kappa系數、交并比(IoU)等指標驗證茶葉分布提取精度。
1)F1值。基于查全率(recall)和查準率(precision)的調和平均,計算公式為:
,
(1)
,
(2)
,
(3)
式中:TP為被模型分類正確的正樣本,即正確分為茶園;FN為被模型分類錯誤的正樣本,即漏分的茶園;FP為被模型分類錯誤的負樣本,即錯分的茶園;TN為被模型分類正確的負樣本,即正確分類的背景樣本。
2)OA計算公式為:
。
(4)
3)Kappa系數。這是一種對遙感圖像的分類精度和誤差矩陣進行評價的多元離散方法,不僅考慮正確分類的像元,同時還考慮了所有錯分、漏分情況,公式為:
,
(5)
。
(6)
4)交并比(IoU),表示A,B2組數據的交集和并集之比,公式為:
,
(7)
式中:A為預測值;B為地面真值。
將不同模型方法獲得的結果以及制作的茶園地面真實標簽進行對比,結果如表1所示,其中白色區域為茶園。從表1中可以看出,本文所提MM-CNN方法可以很好地識別區分出茶園與林地、耕地等其他非茶園地物,充分利用了茶園在Google影像上的紋理特征與時序Sentinel-2影像上茶園光譜特征。采用單一影像的U-Net和1D-CNN方法提取茶園結果中錯分、漏分現象較嚴重。0.5 m分辨率的Google影像上茶園紋理特征明顯,但更新頻率低,并且缺少近紅外波段,對不同綠色植被識別能力較弱,U-Net方法結果中部分耕地(梯田)等其他地物被錯分為茶園。時序Sentinel-2影像光譜特征豐富,可以很好地識別不同類型的綠色植被,但影像空間分辨率相對較低、存在“同物異譜、異物同譜”現象,1D-CNN方法分類結果中碎斑點較多,并且部分林地被錯分為茶園。對比表1中不同方法分類結果可以看出,本文提出的MM-CNN方法分類結果圖斑規則、碎斑點較少、無明顯的錯分現象,但存在部分小面積茶園在U-Net和1D-CNN分類結果被提取出,在MM-CNN方法分類結果中卻被漏分的現象。這主要是Sentinel-2影像空間分辨率(10 m)較低且與Google影像空間分辨率(0.5 m)之間量級差異較大(50倍)導致的,Sentinel-2影像上地物邊緣混合像元效應明顯,對小面積茶園監測能力較弱,二者融合后提取小茶園面積進一步變小,這部分小面積茶園在后處理操作中被識別成碎斑點去除。

表1 不同方法茶園分類結果
為定量評價本文方法的有效性,以野外實調圖斑為參考,利用混淆矩陣計算茶園提取結果的OA和Kappa系數。不同茶園提取方法精度結果如表2所示。由表2可知,本文提出的MM-CNN方法茶園分類結果空間精度高于單一影像源分類結果精度,其茶園提取F1值為0.96,OA為94.10%,Kappa系數為0.86,IoU為0.92。與U-Net結果相比較F1值高出0.12,OA高出14.23%,Kappa系數高出0.27,IoU高出0.20; 與1D-CNN分類結果相比較F1值高出0.09,總體精度高出11.20%,Kappa系數高0.23,IoU高出0.15。

表2 不同茶園提取方法精度結果
本文方法提取的武夷山市新田鎮茶園分布結果如圖9所示,新田鎮茶園在全鎮均有分布,主要集中在新田鎮東北部。根據《武夷山市統計年鑒——2020》[29]中記錄的統計數據,新田鎮2019年茶園面積在2.98萬畝(1)1畝=666.7 m2左右,本文提取面積約2.8萬畝,提取面積精度在90%以上,達到實用要求。提取面積略少于參考面積,這與茶園面積統計方式、研究區地形因素以及可能存在的漏分等多種因素有關。

圖9 研究區茶園分類結果
采用MM-CNN方法提取福建省武夷山市茶園分布,相對于單一影像源的深度學習方法(1D-CNN和U-Net)其優勢體現在:
1)充分利用不同影像源在茶園分布提取中的有效信息,準確提取茶園分布。Sentinel-2影像可見光和近紅外4個波段提供豐富的光譜信息,借助不同地物光譜差異(圖2、圖3)可以識別出茶園。Sentinel-2影像的“異物同譜”現象,使不同地類類間方差減小[14],單期影像茶園分類結果中錯分嚴重。時序Sentinel-2影像光譜差異(圖2、圖3)體現不同植被的物候特征,多期影像組合的方式可以減少“異物同譜”現象帶來的錯分。但茶園種植情況復雜,且茶樹長勢、樹齡不一[14,19],導致茶園在影像上有明顯“同物異譜”現象,使用1D-CNN采用時序Sentinel-2影像提取所有茶園特征需大量樣本,且仍會存在林地與茶園的錯分(表1)。1D-CNN方法適用于樹齡長勢差異小、光譜特征單一的茶園提取。超高分辨率(0.5 m)Google影像具有豐富的紋理信息,U-Net方法通過影像紋理可以區分茶園與林地。研究區內具有條紋狀紋理特征的地物除茶園外,還包括梯田、花卉基地以及竹蓀種植基地等地物,且Google影像缺少近紅外波段對植被識別能力較弱,因此使用U-Net方法采用Google影像的茶園分類結果常與梯田等具條紋狀紋理地物錯分(表1)。同時,由于影像質量等原因,部分茶園在影像上沒有清晰的紋理(圖4(c)),U-Net方法也會將影像上其他地物如作物、部分水域錯分為茶園。U-Net方法適用于規模化種植、打理的茶園提取,且影像上茶園紋理特征與其他地物有明顯的差異。本文提出的MM-CNN方法充分利用了時序Sentinel-2影像地物物候特征以及Google影像紋理信息,采用Google影像茶園與林地不同的紋理信息修正時序Sentinel-2影像林地的錯分,借助時序Sentinel-2影像所提供的不同植被物候信息校正Google影像梯田、水域等地物的錯分。MM-CNN方法通過不同影像結果相互約束的方式準確提取出茶園分布,適用于各種樹齡、長勢情況、種植打理情況的茶園提取。
2)MM-CNN使相對經濟、高效地獲取大范圍高精度茶園分布成為可能。本文在研究區多云多霧、無云質量好遙感影像較難獲取的情況下,采用MM-CNN方法較準確地提取出了茶園分布。基于CNN的遙感影像分類方法其分類精度依賴樣本[30],獲取高精度的分類結果對樣本質量有較高要求,樣本制作費時費力。相對于1D-CNN和U-Net方法,MM-CNN方法對樣本依賴性減弱、魯棒性更強,對單一方法分類結果精度要求較低,僅要求沒有茶園漏分。此外,本文研究區內地物類型豐富,對研究茶園與林地、耕地等背景地物的分類有代表性; 所使用Sentinel-2影像時序較少,Google影像也存在區域色調不一致、鑲嵌痕跡明顯等質量問題,對研究多源影像融合的茶園分類具有普適性。因此MM-CNN方法適用性廣,可為南方丘陵山區大范圍高精度獲取茶園分布提供方法參考。
然而,MM-CNN方法也存在邊緣不準確、小面積茶園會被丟棄的現象。Sentinel-2影像空間分辨率相對較低,影像上地物邊緣混合像元效應明顯,其分類結果邊緣較粗糙,后處理可使邊緣平滑,但邊緣準確性較低。其次Sentinel-2影像空間分辨率為10 m,對小面積茶園監測能力較弱,不易識別出小茶園,Google影像空間分辨率較高為0.5 m,二者分辨率相差50倍,分類結果融合后提取小茶園面積進一步變小,這部分小面積茶園在后處理操作中易被識別成碎斑點去除。除此之外,本研究使用的Google影像由2018—2020年多期影像合成,其時效性較低,可能存在新種植茶園在Google影像上是林地、耕地等茶園開墾前地物特征的情況。因此MM-CNN方法采用時效性低的影像會導致新茶園的漏分。
針對單一影像源茶園難提取問題,本文提出了MM-CNN茶園分類方法,在網絡結構構建上,綜合了1D-CNN與2D-CNN網絡結構。首先結合武夷山市茶園實調數據,分析Sentinel-2影像中對應的光譜信息,構建茶樹重要生長期的物候曲線,應用1D-CNN模型初步提取茶園以及疑似區域(林地); 其次以2D-CNN網絡U-Net為模型,借助不同地物紋理信息對0.5 m分辨率的Google影像中林地和茶園(含耕地等疑似地物)進行提取; 最終融合2個分類結果,得到茶園分布。使用OA和Kappa系數對茶園分類結果精度進行定量評價。實驗結果表明,本文方法茶園分類結果的空間精度高于單一影像分類結果,其茶園提取F1值為0.96,OA為94.10%,Kappa系數為0.86,IoU為0.92。本文提出的MM-CNN方法在進行大面積茶園分布提取時,不僅能有效利用多期影像光譜特征構建的茶樹物候信息區分茶園與其他地物,減少“異物同譜、異物同譜”現象帶來的錯分,而且能有效結合茶園地塊的紋理信息,使最終分類結果精度更高。
此外,盡管本文方法對大范圍茶園種植分布提取有效,但分類結果中存在小面積茶園漏分現象。下一步工作是采用更高分辨率的時序多光譜影像(如GF-2)進行茶園分布提取研究,以期提高小面積茶園的分類精度。