結合光學與雷達影像多特征的熱帶典型人工林樹種精細識別

黃 翀 張晨晨,2 劉慶生 李 賀 楊曉梅 劉高煥

(1.中國科學院地理科學與資源研究所 資源與環境信息系統國家重點實驗室 北京 100101；2.中國科學院大學 北京 100049)

近一個多世紀以來，隨著人類社會經濟的快速發展，大規模天然林采伐已成為全球性問題(Paynetal.，2015)；與此同時，出于對木材的大量需求、經濟利益驅動以及生態保護需要，人工林大面積種植，現已成為陸地生態系統的重要組成部分。東南亞是全球范圍內人工林主要分布和快速擴張地區之一(Paynetal.，2015)，既包括為工業生產提供原材料的速生林種，如桉樹(Eucalyptus)，也包括具有較高經濟價值的經濟作物林種，如橡膠樹(Heveabrasiliensis)、油棕(Elaeisguineensis)。據統計，泰國橡膠林和油棕林面積分別從1961年的約40萬hm2和950 hm2增至2017年的約300萬hm2和76萬hm2，增加近800%和80 000%(http:∥www.fao.org)。人工林建設會對原有的地形地貌和生態環境造成一系列消極影響，如土壤肥力下降、水土流失、生物多樣性降低等(Locatellietal.，2015)，準確獲取人工林類型和空間格局，對森林碳儲量估算、生物多樣性保護以及可持續森林管理規劃等研究至關重要(Georgeetal.，2014)。

遙感在人工林空間分布制圖和時間動態監測方面發揮著重要作用，利用機載高光譜數據(Féretetal.，2012)或激光雷達(LiDAR)數據(Shietal.，2018)能夠獲得高精度的樹種分布信息；然而，由于機載高光譜和激光雷達數據獲取成本高昂，對于較大空間尺度的人工林識別，衛星遙感數據得到了更廣泛應用，許多研究使用光學衛星影像(如MODIS、Landsat、SPOT等)進行人工林監測和制圖，特別是針對橡膠林(Lietal.，2012)、油棕林(Broichetal.，2011)、桉樹林(Somersetal.，2010)等。考慮到人工林冠層尺寸，較低空間分辨率遙感影像對其識別存在一定困難(Wuetal.，2009)，而中高空間分辨率遙感影像的觀測視場與人工林樹種或林分冠層尺寸接近，能夠有效減少混合像元、提高植被的細節表達(Gessneretal.，2013)，研究表明，基于高空間分辨率遙感影像，結合影像光譜特征、紋理特征及其他輔助特征可有效提高樹種分類精度(魏晶昱等，2016；尹凌宇等，2016)。在東南亞熱帶地區，由于雨季長，持續的云層覆蓋嚴重，很難獲得足夠的高質量高空間分辨率光學衛星影像。雷達衛星具有較強穿透性，受云雨天氣影響很小，可全天時全天候工作，在熱帶地區土地覆蓋制圖中得到越來越多重視(Kohetal.，2011)。國內外學者已開始嘗試將光學數據(如MODIS、Landsat)與合成孔徑雷達(synthetic aperture radar，SAR)數據(如PALSAR)相結合進行人工林制圖(Dongetal.，2013；Gutiérrez-Vélezetal.，2013)，根據SAR數據后向散射系數的差異在大類上劃分林地和非林地，人工林與其他林地分類更多依賴光譜特征；但由于這2類衛星傳感器在空間分辨率和影像獲取時相等方面存在較大差異，一定程度上限制了二者協同應用效果。近年來，歐洲新一代哨兵(Sentinel)系列衛星為在高空間分辨率上協同多源數據進行土地覆蓋(變化)監測提供了新途徑，目前，單獨利用Sentinel-1或Sentinel-2數據在土地覆蓋(何云等，2019)、樹種識別(Immitzeretal.，2016)、作物分類(韓濤等，2018)和擾動監測(Hojas-Gasconetal.，2015)等方面已開展了較為深入的探索，但協同Sentinel-1雷達信息和Sentinel-2光學信息的樹種識別研究尚不多見。

鑒于此，本研究以近年來人工林快速擴張的泰國東部地區為例，結合Sentinel-2光譜特征、紋理特征和Sentinel-1雷達后向散射特征，開展光學與雷達影像協同人工林樹種分類研究，評價不同數據源和不同特征在樹種分類中的作用，獲取最優分類策略，以期為熱帶典型人工林樹種精細識別提供新的技術途徑。

1 研究區概況與研究方法

1.1 研究區概況

研究區位于泰國東部Chachoengsao(北柳府)、Chon Buri(春武里府)、Chanthaburi(尖竹汶府)和Rayong(羅勇府)交界處(圖1)，主要土地覆蓋類型包括天然林、人工林、耕地、水體、濕地和不透水面等；屬熱帶季風氣候，日照充足，終年炎熱，年均氣溫28 ℃左右，年均降雨量大于1 000 mm。近年來，由于地方經濟發展需要，大面積天然林或耕地被人工林取代，以橡膠樹、油棕和桉樹為代表的人工經濟林種植面積不斷增加。橡膠樹多種植在平原或坡度較緩的山坡上，行寬株密，行距約8 m，株距約3 m；油棕多種植在平原地區，因種植擴張，近年來向高海拔地區的發展趨勢明顯，采用三角形種植模式，套種間距9 m，以充分利用光照，最大限度提高產量(Basiron，2007)；桉樹人工林排列整齊，大小均勻，樹冠較小，多種植在較為平坦的地區。

1.2 遙感數據與預處理

Sentinel-2光學數據和Sentinel-1雷達數據從歐空局(https:∥scihub.copernicus.eu/)免費下載。Sentinel-2數據獲取時間為2019年1月17日，Sentinel-1數據獲取時間為2019年1月12日。利用歐空局提供的SNAP工具對Sentinel-2和Sentinel-1數據進行預處理。

Sentinel-1數據為Level-1級別地距影像(ground range detected，GRD)IW模式雙極化數據，空間分辨率5 m×20 m，幅寬250 km，包含VH和VV兩種極化模式。對數據進行輻射定標、多視處理、斑點濾波和地形校正，采用Refined Lee散斑濾波器(Zhangetal.，2014)去除斑點噪聲的影響，窗口大小設置為7×7。預處理后Sentinel-1影像的像元亮度轉換為以dB為單位的后向散射系數：

σ0(dB)=10×lg DN。

(1)

式中：σ0為后向散射系數；DN為像元亮度。

Sentinel-2數據覆蓋13個光譜波段，包括10 m空間分辨率的3個可見光波段和1個近紅外波段，20 m空間分辨率的3個紅邊波段、1個近紅外波段和2個短波紅外波段，以及60 m空間分辨率的海岸/氣溶膠、水汽和卷積云波段(本研究未用到)。原始Sentinel-2數據為經輻射校正和幾何校正后的L1C標準產品，利用歐空局提供的Sen2Cor插件進行大氣校正，獲取地表反射率。采用雙線性法將20 m空間分辨率的6個波段重采樣至10 m空間分辨率，最終得到10個10 m空間分辨率的波段。

1.3 Sentinel-2光譜特征計算

為充分利用Sentinel-2的高光譜分辨率，驗證新增紅邊波段在樹種分類中的作用，本研究計算包括15個紅邊指數在內的33個光譜指數(表1)，加上Sentinel-2的10個原始波段，共計43個光譜指數。

1.4 Sentinel-2紋理特征提取

Haralick(1979)提出的灰度共生矩陣(gray-level co-occurrence matrix，GLCM)是一種公認的比較成熟、有效的紋理特征提取方法，其使用一個空間共生矩陣計算像素值之間的關系，并利用這些值計算矩陣的二階統計性質。采用灰度共生矩陣提取紋理特征涉及3個重要參數：窗口大小、步長和移動方向。通過多次試驗對比分析，本研究設置窗口大小為7×7，步長為1個像元，移動方向取0°、45°、90°和135° 4個方向的平均值，提取Sentinel-2紋理特征。選取4個重要的紋理參數，包括角二階矩(angular second moment，ASM)、對比度(contrast，CON)、相關性(correlation，COR)和熵(entropy，ENT)參與樹種分類。

1.5 Sentinel-1雷達后向散射特征提取

不同極化模式的后向散射系數及相應的差值和比值已被證明對影像分類具有重要作用(Dongetal.，2012)。本研究通過下式計算VV和VH的差值(Diff)和比值(Ratio)，選擇VV、VH、Diff和Ratio共4個后向散射特征作為樹種分類的Sentinel-1輸入變量：

(2)

(3)

1.6 特征組合及分類

采用隨機森林(random forest，RF)模型作為樹種分類及特征重要性計算和選擇的工具。RF是Breiman(2001)提出的一種以決策樹為基本分類器的集成學習算法，與其他機器學習算法相比，RF算法能夠并行處理高維、海量數據，對訓練數據中的噪聲或孤立點具有更強的穩定性和魯棒性(Wurmetal.，2017；Rodriguez-Galianoetal.，2012)，且具有較為突出的特征選擇能力，可利用不參與訓練的袋外(out-of-bag，OOB)數據估計每個預測變量的重要性，也可利用OOB數據評估OOB score，以確定RF模型的最佳輸入特征，減少特征冗余。本研究中決策樹數目(ntree)設為800，每個節點的特征數(mtry)由下式計算：

(4)

式中：p為用于分類的特征數。

由于Sentinel-2包含豐富的原始光譜波段及其衍生出的眾多光譜指數特征，全部特征參與分類可能造成信息冗余，導致分類精度降低、分類速度下降，因此，本研究在RF模型中采用平均不純度減少算法(Breiman，2001)對光譜波段和光譜指數進行重要性評估，根據OOB score確定最優光譜特征組合。

為評估Sentinel-2光譜特征、紋理特征和Sentinel-1雷達后向散射特征對人工林樹種識別的能力和貢獻，本研究考慮4種特征組合(表2)作為RF分類器的輸入特征，比較不同類型遙感數據及其特征在分類精度提升上的表現。

表2 分類特征組合Tab.2 Feature combinations for classification

1.7 精度評價

2018年9月，對研究區土地覆蓋類型進行實地調查，獲得251個實地樣本點。通過Google Earth的高分辨率衛星影像(http:∥earth.google.com/)和Global Croplands的實地圖像(https:∥croplands.org)對樣本點進行補充，盡可能使樣本點均勻分布在研究區內。最終共選取1 397個驗證樣本，其中包括203個天然林樣本、207個橡膠林樣本、209個油棕林樣本、131個桉樹林樣本、251個耕地樣本、143個水體樣本、50個濕地樣本和203個建筑用地樣本。采用混淆矩陣對4種特征組合的分類精度進行評價，評價指標包括生產者精度(producer’s accuracy，PA)、用戶精度(user’s accuracy，UA)、總體精度(overall accuracy，OA)、Kappa系數以及PA和UA的調和平均值(F1)。F1的取值范圍為[0，1]，值越大，分類效果越好，反之則越差。F1計算公式(Baumannetal.，2012)如下：

(5)

2 結果與分析

2.1 不同光譜特征在人工林樹種分類中的重要性

采用平均不純度減少算法計算得到Sentinel-2影像15個紅邊指數、18個光譜指數和10個原始波段的特征重要性如圖2所示。由綠波段B3和近紅外波段B8A計算得到的NDWI2的重要性最高，為5.02%；其次為綠波段B3和短波紅外波段B11，重要性分別為4.45%和4.29%；紅邊波段B5和紅邊指數NDVIre1、NDre2和NDVIre2的重要性較高，分別為2.89%、3.97%、3.67%和2.68%，說明紅邊波段對樹種分類具有較高價值；而原始波段B8和B6的重要性最低，不足1%。

圖2 Sentinel-2光譜特征重要性排序Fig.2 Feature importance of spectral features for Sentinel-2 imagery

43個光譜特征組合模型的OOB score分析如圖3所示。特征數從1增至17，OOB score逐漸增加，特征數為17時達最高值0.947 2，此后有輕微降低趨勢，因此本研究選擇重要性排名前17的光譜特征參與樹種分類。統計重要性排名前17的光譜特征中包含Sentinel-2各波段的特征數，如圖4所示。在前17個光譜特征中，包含近紅外波段B8A的特征最多(8個)，其次為包含可見光波段B3(7個)、B4(6個)的特征。此外，包含紅邊波段B5的特征也較多，進一步說明新增紅邊波段在樹種分類中具有較高價值。

圖3 Sentinel-2不同光譜特征組合模型的OOB scoreFig.3 OOB score of different feature combination models for Sentinel-2 imagery

圖4 重要特征波段統計Fig.4 Statistics of selection for each Sentinal-2 band

2.2 不同人工林樹種在紋理特征和后向散射特征上的可分離性

分別提取天然林、桉樹林、油棕林和橡膠林在紋理特征和后向散射特征上的特征值(圖5)，4種林地均具有一定可分離性。從Sentinel-2紋理特征(圖5a)看，天然林在ASM特征上明顯小于人工林，在CON和ENT特征上明顯大于人工林，可分離性較好；在COR特征上與人工林區分度較小。對于3種人工林，油棕林在CON和ENT特征上值最大，在ASM特征上值最小，COR特征值介于桉樹林和橡膠林之間；桉樹林在ASM、COR特征上值最大，在ENT特征上值最小；橡膠林在CON和COR特征上值最小，ASM和ENT特征值介于桉樹林和油棕林之間。這說明桉樹林、油棕林和橡膠林在4種紋理特征上具有一定差異性。從Sentinel-1雷達后向散射特征(圖5b)看，4種林地均具有一定可分離性，在Ratio特征上的區分度大于Diff特征。總的來說，Sentinel-2紋理特征和Sentinel-1雷達后向散射特征可作為樹種分類的有效特征。

圖5 天然林和不同人工林在紋理特征(a)和后向散射特征(b)上的特征值(均值±標準差)Fig.5 Statistics of textural(a)and backscattering features(b)of natural forest and different plantations(mean±SD)圖5a的主坐標軸表示ASM、CON和COR特征值，次坐標軸表示ENT特征值；圖5b的主坐標軸表示VV、VH和Diff特征值，次坐標軸表示Ratio特征值。The primary axis of Fig.5a represents the values of ASM,CON and COR features,and the secondary axis represents the values of ENT feature;the primary axis of Fig.5b represents the values of VV,VH and Diff features,and the secondary axis represents the values of Ratio feature.

2.3 不同特征組合樹種分類結果比較

4種特征組合的分類結果如圖6所示。僅利用Sentinel-2光譜特征(S)分類(圖6a)，“椒鹽現象”非常明顯，加入紋理特征(S+T)、后向散射特征(S+SAR)以及結合紋理特征和后向散射特征(S+T+SAR)(圖6b-d)后，各地類分類圖斑的破碎度均有所降低，“椒鹽現象”得到改善。3種人工林的分類圖斑較為規整，其中，橡膠林分布最連續，面積最大，主要分布在平原或坡度較緩的山坡上；油棕林多與橡膠林鑲嵌分布，地塊面積較大；桉樹林地塊面積最小，主要分布在較為平坦的地區，多與耕地交錯分布。

圖6 不同特征組合的分類結果Fig.6 Classification results of different feature combinations

不同特征組合分類結果的混淆矩陣如圖7所示，精度驗證結果如圖8所示。對比4種特征組合的分類精度，僅利用光譜特征(S)分類，3種人工林的分類精度均較低，不同人工林之間存在不同程度混淆，且易與天然林和耕地混淆；橡膠林和油棕林的生產者精度不足0.70，桉樹林的生產者精度僅0.53，F1僅0.61，總體分類精度為0.75，Kappa系數為0.71。加入紋理特征組合(S+T)分類，天然林與3種人工林之間混淆減少，人工林誤分為耕地比例減小，桉樹林、油棕林和橡膠林的生產者精度分別提高0.16、0.15和0.08；3種人工林的整體分類精度均有所提升，桉樹林、油棕林和橡膠林的F1分別提高0.16、0.11和0.04，總體分類精度提高至0.80，Kappa系數提高至0.77。加入后向散射特征組合(S+SAR)分類，桉樹林、油棕林和橡膠林相比僅利用光譜特征(S)的分類精度同樣有所提升，F1分別提高0.18、0.15和0.08，總體分類精度提高至0.84，Kappa系數提高至0.81。結合Sentinel-2光譜特征、紋理特征和Sentinel-1雷達后向散射特征(S+T+SAR)分類，總體分類精度和Kappa系數達到最高值，分別為0.85和0.83；其他地類誤分為人工林的數量最少，3種人工林均達到最高用戶精度；桉樹林、油棕林和橡膠林的F1均大于0.80，相比僅利用光譜特征(S)分類分別提高0.19、0.15和0.10。

圖7 不同特征組合分類結果的混淆矩陣Fig.7 Confusion matrix of classification results of different feature combinations

圖8 不同特征組合分類精度對比Fig.8 Accuracies of classification results of different feature combinations

綜上可知，4種特征組合的分類結果中，S+T+SAR的總體分類精度和Kappa系數最高，其次為S+SAR和S+T，S最低。桉樹林、油棕林、橡膠林、耕地、水體和建筑用地在S+T+SAR中均取得最好分類結果，天然林和濕地在S+SAR中分類精度最高。

3 討論

本研究基于Sentinel-2光學影像豐富的光譜特征、紋理特征并結合Sentinel-1雷達后向散射特征，開展熱帶典型人工林樹種精細識別，結果發現，Sentinel-2的短波紅外波段B11在樹種分類中重要性較高，與Immitzer等(2016)和Sothe等(2017)基于Sentinel-2A進行樹種分類時得到的結論一致。在隨機森林算法選擇的17個光譜特征中，Sentinel-2的紅邊波段和紅邊指數所占比例較大，許多學者(Immitzeretal.，2016；Sotheetal.，2017；Mngadietal.，2019)同樣強調了紅邊波段對林地分類的重要性。這是由于植被在紅邊波段敏感度高(Peerbhayetal.，2014)，尤其是紅邊波段與葉綠素含量、生物量和冠層結構等反映植被葉片性質的信息密切相關(Dubeetal.，2014；Ramoeloetal.，2015)，因此對不同樹種的區分具有重要意義。此外，Sentinel-2光學影像的紅波段和近紅外波段在樹種分類中同樣表現出較高重要性，這可能是因為紅波段葉綠素吸收能力強，近紅外波段吸水能力強，從而提高了不同林地之間及與周圍地物的光譜可分性。Nooni等(2014)基于Landsat-7提取油棕林時也發現在紅波段和近紅外波段油棕的可分離性最高。

僅利用Sentinel-2光譜特征分類，橡膠林、油棕林和桉樹林3種人工林的識別精度均較低。由于光譜特征相似，人工林與天然林、耕地混淆較多，不同人工林之間也存在不同程度混淆。Dian等(2015)研究表明，結合光譜特性和紋理特征可有效解決光譜混淆造成的分類誤差，提高樹種分類精度。本研究在光譜特征基礎上加入紋理特征(S+T)后，3種人工林的分類精度均有不同程度提高。不同人工林具有不同的株行距和不同形狀、大小的樹冠，在遙感影像上會形成特有的紋理，而天然林在樹冠大小、樹高和密度等方面具有更大變異性，紋理特征加入可有效減少“椒鹽現象”，提高斑塊的完整性和分類精度。如僅利用光譜特征分類，緩坡上的橡膠林因海拔、坡向、坡度等因素影響呈現不同的光譜特征，易被錯分為耕地和桉樹林；加入紋理特征后，可有效降低山體效應帶來的光譜差異，減少橡膠林與桉樹林間的誤分，橡膠林斑塊也更為完整。

加入后向散射特征(S+SAR)后，3種人工林的分類精度同樣有所提升，說明Sentinel-1雷達后向散射特征可有效提高人工林樹種間的區分度。Torbick等(2016)在提取緬甸和加里曼丹的橡膠林和油棕林時發現，Senitnel-1雷達后向散射特征對橡膠林和油棕林提取具有重要作用。雷達數據基于地物后向散射特征獲得不同于光學遙感的影像，且雷達信號具有一定穿透力，可獲取植被表面信息和地表粗糙度等額外信息，尤其是合成孔徑雷達(SAR)對橡膠林、油棕林和桉樹林在生物量、密度和垂直分層等森林結構信息上的差異敏感性較高，使得SAR數據在樹種識別中具有獨特優勢。此外，不同樹種冠層含水量也可通過影響介電常數影響雷達信號(Oonetal.，2018)，能夠提高不同人工林樹種間的區分度。因此，SAR數據可作為光學影像的有益補充，對光譜特征接近的人工林和作物進行區分。研究區耕地多種植水稻(Oryzasativa)，淹水期長，SAR后向散射特征對土壤濕度和淹沒程度十分敏感(Chatziantoniouetal.，2017)，可以彌補光學遙感區分人工林和水稻的不足，提高人工林分類精度。

近年來研究發現，集成多種觀測類型和數據模式的多傳感器與數據融合技術有助于提高森林制圖的細節和精度(Hydeetal.，2006；Descleeetal.，2013)。本研究中，僅利用Sentinel-2光譜特征分類精度較低，加入光學紋理特征和雷達后向散射特征(S+T+SAR)后，3種人工林的分類精度均達到最高，這說明單一傳感器或特征進行人工林樹種分類仍有很大局限，尤其在生態系統復雜、植被覆蓋度高的東南亞地區；而多數據源、多特征結合可充分利用地物在不同傳感器上感知的信息，增加不同人工林之間以及人工林與其他地物之間的區分度。本研究結合Sentinel-2光學數據和Sentinel-1雷達數據的人工林樹種分類精度明顯優于以往僅利用Sentinel-2數據的樹種分類精度(Immitzeretal.，2016)，可見多傳感器、多特征結合應用的重要性。但需要說明的是，本研究Sentinel-1數據只提取了后向散射特征，實際上由于雷達對植被冠層的穿透能力和對土壤含水量的敏感性，高空間分辨率Sentinel-1數據同樣含有豐富的紋理信息。此外，Sentinel-1、Sentinel-2數據均具有較高的重訪周期，如何結合其多樣的紋理特征和時序特征開展精細化人工林制圖值得進一步研究。

4 結論

本研究以近年來人工林快速擴張的泰國東部地區為例，結合Sentinel-2光譜特征、紋理特征和Sentinel-1雷達后向散射特征，采用隨機森林算法對Sentinel-2不同光譜特征在樹種分類中的重要性進行評估，并分別針對光譜特征、紋理特征和后向散射特征的不同組合進行3種人工林樹種識別研究，最終實現該地區3種人工林樹種的精細提取，得到如下結論：

1)在Sentinel-2光學影像的43個光譜特征中，選擇重要性排名前17的特征參與分類時，OOB score達最高值0.947 2；Sentinel-2的藍波段、紅邊和近紅外波段及其相應的植被指數在樹種識別中重要性較高。

2)桉樹林、油棕林、橡膠林和天然林在Sentinel-2影像的ASM、CON、COR和ENT紋理特征和Sentinel-1影像的VV、VH、Diff和Ratio后向散射特征上均具有一定差異性，可作為樹種分類的有效特征。

3)僅利用Sentinel-2光譜特征對不同人工林樹種的區分度有限，桉樹林、油棕林和橡膠林的F1分別為0.61、0.74和0.70。結合光譜特征、紋理特征和后向散射特征，3種人工林的分類精度達到最高，F1均大于0.80，同時其他地物的分類精度也有較大提高，總體分類精度和Kappa系數均達到最高值，分別為0.85和0.83。