劉靈,張加龍*,韓雪蓮,許冬凡,王書賢,程滔
(1.西南林業(yè)大學(xué) 林學(xué)院,昆明 650224;2.國家基礎(chǔ)地理信息中心,北京 100830)
作為森林資源調(diào)查的重要內(nèi)容,森林樹種的識別與分布情況對于森林資源管理、生物多樣性評估、水土保持以及森林碳儲量估算都具有重要參考價值[1]。
與傳統(tǒng)的野外森林調(diào)查相比,衛(wèi)星遙感技術(shù)因其快速和高效被廣泛應(yīng)用于森林分類與識別中。各類型的遙感影像被用于樹種識別[2-4]。有研究表明,數(shù)據(jù)源也是影響樹種識別精度的因素之一[5]。多數(shù)情況下,中低空間分辨率主要用于森林類型分類,而未用于樹種組成分析;高光譜數(shù)據(jù)和激光雷達(dá)數(shù)據(jù)因其高昂的成本而無法應(yīng)用于大尺度區(qū)域樹種識別[6-7]。Sentinel-1雷達(dá)數(shù)據(jù)和Sentinel-2多光譜數(shù)據(jù)憑借其免費獲取和較高空間分辨率的優(yōu)勢在大區(qū)域樹種精細(xì)識別領(lǐng)域被廣泛使用。
雷達(dá)和光學(xué)影像的融合可以克服單個傳感器的缺點,充分利用多種信息之間的優(yōu)勢互補來提高地物提取精度,尤其改善森林分類[8]。黃翀等[8]通過Sentinel-1/2數(shù)據(jù)提取樹種的多種特征,利用隨機森林算法進(jìn)行熱帶典型人工林樹種精細(xì)識別,總體精度達(dá)85%。有研究表明,多時相時間序列數(shù)據(jù)能反映不同樹種在生長發(fā)育過程中隨時間的周期性變化,識別不易區(qū)分的樹種類型,提高樹種分類的精度[9]。例如,Bjerreskov等[10]基于多時相Sentinel-1和Sentinel-2數(shù)據(jù)結(jié)合隨機森林算法,對丹麥景觀進(jìn)行全面分類,結(jié)果表明多時相和多傳感器數(shù)據(jù)的結(jié)合產(chǎn)生了較好的分類效果。Ghorbanian等[11]通過構(gòu)建Sentinel-1和Sentinel-2時間序列數(shù)據(jù)以提取樹種的季節(jié)物候特征,并利用隨機森林算法進(jìn)行紅樹林生態(tài)系統(tǒng)制圖,總體精度達(dá)93.23%,研究結(jié)果表明季節(jié)特征和多源遙感數(shù)據(jù)的結(jié)合有助于獲得更高的分類精度。
除數(shù)據(jù)源外,分類算法在一定程度上也會影響分類結(jié)果。隨機森林分類算法因其模型訓(xùn)練少,計算精度高而被廣泛應(yīng)用于森林類型分類中,是目前森林樹種識別應(yīng)用最廣泛的分類方法[8,10-11]。
進(jìn)行大范圍森林樹種識別時,數(shù)據(jù)的獲取與處理仍是關(guān)鍵問題。Google Earth Engine(GEE)平臺憑借其并行運算,不受空間和時間約束的優(yōu)勢,在土地利用分類[12]、植被類型分類[13]和作物識別[14]等方面被用于進(jìn)行深入研究。同時,研究區(qū)森林類型以針葉林為主,占全市面積的80%以上,在香格里拉市生態(tài)保護以及社會經(jīng)濟發(fā)展中占據(jù)重要地位[15]。因此,本研究基于GEE平臺,以Sentinel-1/2時間序列數(shù)據(jù)為主要數(shù)據(jù)源,充分利用光學(xué)數(shù)據(jù)的光譜和紋理特征以及合成孔徑雷達(dá)(SAR)數(shù)據(jù)的雷達(dá)和紋理特征,結(jié)合時序特征和地形特征,使用隨機森林算法對香格里拉4種針葉林優(yōu)勢樹種識別進(jìn)行研究,以期為該區(qū)域進(jìn)一步開展森林碳儲量、碳匯、生態(tài)水文保護等研究提供科學(xué)參考。
香格里拉市位于云南省西北部的滇川藏三省交界區(qū),地處迪慶香格里拉腹心地帶,地理坐標(biāo)為99°20′~100°19′E,26°52′~28°52′N,如圖1所示。研究區(qū)地形西北高、東南低,平均海拔3 549 m,垂直落差大,氣候?qū)贉貛?山地季風(fēng)氣候,年平均氣溫4.7~16.5 ℃,是我國保存較好的以亞高山針葉林為主的天然林區(qū)。據(jù)統(tǒng)計,研究區(qū)總面積11 613 km2,森林覆蓋率74.99%,植物資源豐富,常見樹種有高山松(Pinusdensata)、云南松(Pinusyunnanensis)、云冷杉(Piceaasperata)和高山櫟(Quercusaquifolioides)等[15]。
1.2.1 Sentinel數(shù)據(jù)
Sentinel-1雷達(dá)數(shù)據(jù)和Sentinel-2光學(xué)數(shù)據(jù)均從GEE平臺獲取(https://developers.google.com/earth-engine/datasets/catalog/sentinel),獲取的影像時間范圍均為2020年1月1日—2020年12月31日。
Sentinel-1數(shù)據(jù)為IW模式下VV、VH雙極化特征的Level-1級別地距影像(Ground Range Detected,GRD),分辨率為10 m。GEE平臺中的Sentinel-1 SAR GRD數(shù)據(jù)已經(jīng)經(jīng)過Sentinel-1工具箱工具進(jìn)行了多視處理、應(yīng)用軌道文件和輻射定標(biāo)等預(yù)處理。由于SAR影像中固有的斑點噪聲降低了影像的質(zhì)量,因此本研究采用Refined Lee散斑濾波器消除斑點噪聲的影響[16]。
Sentinel-2數(shù)據(jù)以高空間分辨率和豐富的光譜信息被廣泛應(yīng)用于植被、土壤和海洋等監(jiān)測和分類領(lǐng)域。本研究選擇已經(jīng)過大氣校正的Sentinel-2A數(shù)據(jù),選取B2—B8、B8A、B11、B12共計10個波段,并用雙線性法將波段重采樣至10 m。由于云及陰影的噪聲影響,篩選云量低于30%的影像中值合成無云影像。
圖1 研究區(qū)位置、樣本及影像Fig.1 Study area, samples and images
主成分分析(Principal Component Analysis,PCA)是一種用于消除數(shù)據(jù)冗余并保留相關(guān)信息的技術(shù)。使用PCA對SAR與光學(xué)數(shù)據(jù)進(jìn)行融合不會丟失大部分信息,因此,與其他融合方法相比,基于PCA融合能提供更高質(zhì)量的數(shù)據(jù)[13]。本研究首先對Sentinel-2影像進(jìn)行PCA變換,然后將貢獻(xiàn)最少信息量的2個成分(最后2個成分)替換為Sentinel-1的VV和VH波段,最后通過逆主成分變換將波段融合,得到數(shù)據(jù)質(zhì)量較高的影像。
1.2.2 DEM(數(shù)字高程模型)數(shù)據(jù)
研究區(qū)地形垂直落差大,樹種的分布具有垂直分布的特點,結(jié)合地形數(shù)據(jù)能夠有效增強樹種的可分性[17-18]。本研究所使用的高程數(shù)據(jù)(SRTM DEM)直接從GEE獲取,空間分辨率為30 m,并生成海拔(Elevation)、坡度(Slope)、坡向(Aspect)。
1.2.3 樣本數(shù)據(jù)
以Google Earth高清歷史影像和森林資源調(diào)查數(shù)據(jù)為參考繪制樣本數(shù)據(jù)。森林資源調(diào)查數(shù)據(jù)來源于香格里拉市2016年森林資源二類調(diào)查成果,以矢量形式存儲,包括胸徑、樹高、地類、優(yōu)勢樹種和面積等屬性信息。
本研究使用的樣本數(shù)據(jù)共分為3層。第1層樣本數(shù)據(jù)為通過Google Earth高清影像目視解譯獲得的300 m×300 m面數(shù)據(jù),用于森林和非森林分類;第2層樣本數(shù)據(jù)由森林資源調(diào)查數(shù)據(jù)獲得,根據(jù)地類和優(yōu)勢樹種隨機生成30 m×30 m的面數(shù)據(jù),用來劃分針葉林和闊葉林。第3層樣本數(shù)據(jù)來源于森林資源調(diào)查數(shù)據(jù),根據(jù)優(yōu)勢樹種隨機生成10 m×10 m的面數(shù)據(jù),用來劃分優(yōu)勢樹種。每層數(shù)據(jù)均由面數(shù)據(jù)轉(zhuǎn)為點數(shù)據(jù),并選取70% 作為訓(xùn)練數(shù)據(jù),30% 作為測試數(shù)據(jù),訓(xùn)練數(shù)據(jù)和測試數(shù)據(jù)各自獨立不重復(fù),樣本統(tǒng)計見表1。
表1 分類體系及樣本數(shù)量
優(yōu)勢樹種識別技術(shù)路線圖如圖2所示。研究過程總體分為3部分,即數(shù)據(jù)獲取與預(yù)處理、特征提取與特征組合以及分層分類與精度評價。首先,在GEE平臺上獲取Sentinel-1/2數(shù)據(jù)并進(jìn)行去噪、去云、裁剪和融合等預(yù)處理;其次,提取Sentinel-1雷達(dá)特征、紋理特征,Sentinel-2光譜特征、紋理特征、時序特征和地形特征等共計43個特征,根據(jù)不同特征組合形成6種組合方案;最后,利用隨機森林算法在3個層次上對6種組合方案分別進(jìn)行森林分類,并對分類結(jié)果進(jìn)行精度評價。
圖2 優(yōu)勢樹種識別技術(shù)路線圖Fig.2 Flowchart of dominant species classification
2.1.1 Sentinel-1雷達(dá)特征提取
研究表明,VV、VH極化特征的差值(Diff)和比值(Ratio)對分類具有重要作用[19]。選擇VV、VH、Diff、Ratio4個波段作為Sentinel-1雷達(dá)特征參與分類,Diff和Ratio計算公式如下
(1)
(2)
2.1.2 Sentinel-1紋理特征提取
由于雷達(dá)對土壤含水量敏感且對植被冠層具有很強的穿透能力,因此Sentinel-1雷達(dá)數(shù)據(jù)包含較豐富的紋理信息。灰度共生矩陣(Gray-Level Co-occurrence Matrix,GLCM)是Haralick[20]提出的一種有效提取紋理特征的方法,其基于像素點對間的灰度關(guān)系統(tǒng)計紋理信息。為減少特征冗余,選取6個重要的紋理特征參與分類,包括熵(Entropy,ENT)、相關(guān)性(Correlation,COR)、對比度(Contrast,CON)、角二階矩(Angular Second Moment,ASM)、相異性(Dissimilarity,DISS)和方差(Variance,VAR)。
2.1.3 Sentinel-2光譜特征提取
研究使用的光譜波段分別為B2—B8、B8A、B11、B12。同時,植被指數(shù)能較好地反映植被生長狀況,因此加入歸一化植被指數(shù)(Normalized Difference Vegetation Index,NDVI)和增強型植被指數(shù)(Enhanced Vegetation Index,EVI)。
2.1.4 時序特征提取
現(xiàn)有的研究表明,植被指數(shù)時間序列被廣泛地應(yīng)用于分類和時序特征提取中[21]。因此,本研究通過重構(gòu)NDVI、EVI時間序列提取植被的時序特征參與分類。本研究用到的時間序列重構(gòu)方法為HANTS諧波分析。其核心原理是時間波譜數(shù)據(jù)可以分解為多個簡單的波譜曲線,首先通過快速傅里葉變換得到觀測樣本的原本擬合曲線,之后通過對無效數(shù)據(jù)賦予零權(quán)重而去除無效數(shù)據(jù),從而將原來的曲線擬合問題轉(zhuǎn)化為加權(quán)最小二乘曲線擬合以達(dá)到時序重構(gòu)的目的。表達(dá)式如下
(3)
式中:A0為余項,表示時序的平均值;n為時間間隔;N為時間序列長度;An為振幅;ωn為時間頻率;φn為相位角。
與AG擬合、DL擬合、Savitzky-Golay濾波相比,HANTS諧波分析法能更好地反映植被在長時間序列內(nèi)的物候特征[20]。其中,相位角(Phase)、振幅(Amplitude)和頻率(Frequency)是HANTS諧波分析重構(gòu)時間序列曲線的主要特征參數(shù)。根據(jù)植被變化規(guī)律,本研究選取的時序特征為相位角(NDVI_ Phase、EVI_ Phase)、振幅(NDVI_ Amplitude、EVI_ Amplitude)和頻率(NDVI_ Frequency、EVI_ Frequency)。
2.1.5 Sentinel-2紋理特征提取
同樣地,Sentinel-2含有豐富的紋理特征。為了避免特征冗余,本研究選用3×3窗口,利用PCA進(jìn)行數(shù)據(jù)降維,基于第一主成分提取6種紋理特征參與分類,詳見表2。
2.1.6 特征組合
本研究設(shè)計了6種特征組合方案來評估不同特征對香格里拉市優(yōu)勢樹種識別的能力與貢獻(xiàn),見表3。
表2 特征統(tǒng)計與說明
表3 特征組合
Breiman[22]提出的隨機森林(Random Forest,RF)分類模型是一種常用的機器學(xué)習(xí)分類算法,其通過構(gòu)建多個決策樹來運行。與其他機器學(xué)習(xí)算法相比,RF算法具有計算精度高、模型訓(xùn)練時間少,且能夠確定變量在模型中的重要性等優(yōu)勢,因此該方法被廣泛應(yīng)用于分類研究中。在GEE中調(diào)用隨機森林分類模型,該模型參數(shù)較少,決策樹數(shù)量和變量數(shù)量是影響其結(jié)果的參數(shù)。針對每一層分類,在計算可以負(fù)荷的情況下,選擇精度最高時的決策樹數(shù)量。研究采用的隨機森林決策樹數(shù)量第1層是100,第2層是300,第3層是200,且變量數(shù)量默認(rèn)為特征變量總數(shù)的平方根。
本研究使用混淆矩陣對分類精度進(jìn)行評價,即用戶精度(User Accuracy,UA,公式中為UA)、制圖精度(Producer Accuracy,PA,公式中為PA)、總體精度(Overall Accuracy,OA)和Kappa系數(shù)。但在實際分類中,PA和UA又相互制衡,不能準(zhǔn)確判斷分類的優(yōu)劣,所以引入PA和UA的調(diào)和平均值(F1)評價分類效果[23]。F1的計算公式如下
(4)
式中,F(xiàn)1的取值范圍是[0,1]。
基于不同特征組合的森林/非森林分類精度評價結(jié)果見表4。所有實驗的OA均高于87%,Kappa系數(shù)均高于0.80,F(xiàn)1均高于86%?;诜桨?所有特征的分類精度最高,OA為98.73%,Kappa系數(shù)為0.97,F(xiàn)1為98.71%。采用精度最高的方案6的分類結(jié)果計算林地面積,得到分類結(jié)果:森林面積為7 856.70 km2,而森林資源調(diào)查面積為7 800.26 km2,兩者之間相對誤差0.72%。通過目視評價和面積對比分析可知,基于方案6的森林/非森林分類結(jié)果可以作為第2層分類的掩膜。
表4 森林/非森林分類精度評價Tab.4 Accuracy assessment for the forest/non-forest classification
基于不同特征組合的針葉林/闊葉林分類精度評價結(jié)果見表5。所有實驗的OA均高于83%,Kappa 系數(shù)均高于0.70,F(xiàn)1均高于83%?;诜桨?所有特征的分類精度最高,OA為92.80%,Kappa系數(shù)為0.85,F(xiàn)1為92.58%。采用精度最高的方案6的分類結(jié)果計算針葉林面積,得出研究區(qū)針葉林面積6 437.94 km2,而森林資源調(diào)查為6 157.42 km2,兩者之間相對誤差4.56%,此次針葉林/闊葉林分類結(jié)果可以作為第3層分類的掩膜。
表5 針葉林/闊葉林分類精度評價Tab.5 The summary of the accuracy assessment for the needle-leaved forest/broad-leaved forest classification
分別提取針葉林中高山松、云南松、云冷杉和落葉松在光譜特征和雷達(dá)特征上的特征均值,如圖3所示。由圖3(a)可知,同為針葉林,不同樹種在不同波段的反射率不同,但在整體上這4種樹種的曲線走勢較為相似,符合綠色植被的光譜曲線。在B2、B3、B4波段呈現(xiàn)“藍(lán)谷”“綠峰”“紅谷”現(xiàn)象,之后波段反射率先上升后下降,在B5、B6、B7、B8、B8A、B11、B12波段相差較大,從而提高不同樹種之間的光譜可分性。由圖3(b)可知,高山松在VV、VH波段與其他3種樹種具有明顯差異,特征值最低,可分離性較好;云南松、云冷杉、落葉松在VV、VH、Diff特征上有一定的差異,且特征值從大到小排序為:云南松、云冷杉、落葉松,可充分利用Sentinel-1雷達(dá)波段對優(yōu)勢樹種進(jìn)行有效提取。
NDVI、EVI時間序列曲線如圖3(c)和圖3(d)所示。云南松的EVI值要高于其他3種類型,NDVI值大部分時間高于其他3種類型。云南松與云冷杉在NDVI時間序列曲線上呈現(xiàn)季節(jié)性相反狀態(tài)。落葉松屬落葉樹種,在冬季的NDVI和EVI值明顯低于其他3種樹種,春夏2季樹木生長,其NDVI和EVI值上升,秋季開始落葉,NDVI和EVI值開始下降。4種樹種的EVI整體呈現(xiàn)春夏上升,秋冬下降,夏季值最高的趨勢。4種樹種的植被指數(shù)時間序列具有一定的差異,可以作為樹種識別的分類基礎(chǔ)。
圖3 4種樹種在光譜特征、雷達(dá)特征和時序特征上的可分離性Fig.3 Separability of spectral features, SAR features, time-series features of 4 species
基于6種特征組合的樹種分類結(jié)果如圖4所示。與僅利用Sentinel-2光譜特征(S)分類相比,加入其他分類特征(S+TS、S+T2、S+SAR、S+SAR+T1、S+TS+T2+SAR+T1+D)后,分類圖斑的破碎化程度減少,各地類更為集中。但是不同特征組合分類結(jié)果的局部差異仍然存在,方案6(如圖5(f)所示)的分類結(jié)果顯示各樹種更聚集,高山松和云南松之間、云冷杉和落葉松之間的混淆程度相對其他方案(圖5(a)—圖5(e))減少,與實際調(diào)查情況(圖5(g))最接近。
圖4 不同特征組合的樹種分類結(jié)果Fig.4 Classification results of different feature combinations
結(jié)合不同特征組合的混淆矩陣(圖6)和精度評價(表6)可知,所有實驗的OA均高于81%,Kappa系數(shù)均高于0.75,F(xiàn)1均高于80%。方案1中,不同樹種之間存在不同程度的混淆,整體分類效果不佳,這是因為各樹種同屬森林地類,僅利用光譜信息很難進(jìn)行區(qū)分,OA為81.28%,Kappa系數(shù)為0.75,F(xiàn)1為80.74%。方案2—方案6中,在光譜特征的基礎(chǔ)上,加入其他分類特征,使得高山松和云南松、云冷杉和落葉松之間的混淆程度有所減少,各優(yōu)勢樹種的PA和UA明顯提高,與方案1相比OA分別提高了2.05%、1.23%、0.82%、2.05%和8.23%,Kappa系數(shù)分別提高了2.68%、1.65%、1.00%、2.72%和11.13%,F(xiàn)1分別提高了1.75%、1.34%、0.75%、2.16%和8.62%,說明增加時序特征、紋理特征、雷達(dá)特征和地形特征在一定程度上可提高樹種識別精度。在方案6中,4種優(yōu)勢樹種間的混淆程度最小,分類結(jié)果最好,整體分類精度最高,其OA為89.51%,Kappa系數(shù)為0.86,F(xiàn)1為89.36%。
綜合來看,高山松和云南松之間容易發(fā)生混淆,原因是:雖然高山松在垂直分布上較云南松高,但在海拔3 000 m以下以高山松與云南松混交為主,因此容易相互混淆。云冷杉和落葉松相互混淆也較多,這是因為這2種優(yōu)勢樹種同屬于寒溫性針葉林,在海拔上分布幾乎一致,而落葉松往往是云冷杉經(jīng)火燒采伐后恢復(fù)更新的樹種,因此比較容易造成錯分。
(A)、(B) 高清影像;(a) S;(b) S+TS;(c) S+T2;(d) S+SAR(e) S+SAR+T1;(f) S+TS+T2+SAR+T1+D;(g) 二類調(diào)查
Pd代表高山松;Py代表云南松;Pa代表云冷杉;Lg代表落葉松。
表6 不同特征組合的分類精度Tab.6 Classification accuracy of different feature combinations %
基于方案6使用RF計算得出的特征重要性結(jié)果見表7。
表7 基于方案6特征重要性結(jié)果Tab.7 The result of feature importance based on experiment 6
海拔的重要性最高,為6.10%;其次為短波紅外波段B11、B12,重要性分別為3.98%和3.33%;紅邊波段B8A、B5、B7、B6的重要性分別為2.64%、2.51%、2.46%和2.32%,近紅外波段B8的重要性為2.45%,也進(jìn)一步說明了Sentinel-2新增紅邊波段、近紅外波段與短紅外波段在樹種識別中的優(yōu)勢。植被指數(shù)EVI和NDVI的重要性為2.73%和2.63%,時序特征EVI_Phase、NDVI_Phase、NDVI_Frequency、NDVI_Amplitude、EVI_Frequency和EVI_ Amplitude的重要性分別為3.19%、2.73%、2.62%、2.47%、2.40%和2.17%,說明了加入植被指數(shù)及其衍生的時序特征能夠有效提升樹種識別的能力;雷達(dá)特征Diff、Ratio、VV和VH的重要性分別為3.13%、2.33%、2.26%和2.17%,說明雷達(dá)特征對樹種識別具有較高價值;Sentinel-2紋理特征pc1_CON、pc1_VAR、pc1_DISS和pc1_COR的重要性分別為2.81%、2.50%、2.28%和2.27%,而波段pc1_ENT和pc1_ASM的重要性最低,不足1%。
本研究基于GEE平臺,以香格里拉市為研究區(qū),通過融合Sentinel-1和Sentinel-2時間序列數(shù)據(jù),提取其時序特征,并結(jié)合雷達(dá)特征、紋理特征和地形特征,采用隨機森林算法和分層分類法開展研究區(qū)優(yōu)勢樹種識別研究,主要得出以下結(jié)論。
(1)融合Sentinel-1/2時間序列數(shù)據(jù)能夠在10 m空間分辨率下進(jìn)行森林樹種識別,總體精度可達(dá)89.51%,Kappa系數(shù)可達(dá)0.86,表明Sentinel數(shù)據(jù)在識別植被這一領(lǐng)域有較好的應(yīng)用價值。
(2)不同分類特征的結(jié)合比僅用光譜特征分類時的精度都高,說明單一特征進(jìn)行樹種識別仍具有局限性,而多源遙感數(shù)據(jù)融合與多特征的結(jié)合能充分利用不同樹種在各個特征上的可分性,提高其識別精度。
(3)分層分類法和隨機森林算法的結(jié)合能夠提高分類效率,使結(jié)果更加可靠。
(4) 研究區(qū)氣候類型多變,植被復(fù)雜多樣,云與地形等因素嚴(yán)重影響遙感影像的獲取,所以研究使用密集的時序數(shù)據(jù)結(jié)合HANTS時序重構(gòu)方法以減少云對識別精度的影響。時間序列影像數(shù)據(jù)及其時序特征對于提高識別精度具有關(guān)鍵作用。
本研究融合Sentinel-1/2時間序列影像的森林優(yōu)勢樹種識別精度明顯高于使用單期Sentinel-1/2影像的樹種分類精度[8],表明采用多源時序遙感數(shù)據(jù)的時序特征能夠提高森林優(yōu)勢樹種的分類精度。與其他使用時序特征(平均值)的樹種識別研究相比[10],本研究加入了基于HANTS時序諧波重構(gòu)的時序擬合特征(相位角、振幅和頻率),分類精度由74%提高到89%。但需說明的是,本研究僅采用了2020年的NDVI、EVI曲線衍生的時序特征進(jìn)行分類,實際上Sentinel-2包含豐富的光譜信息及衍生的光譜指數(shù),如何利用多樣的長時序特征結(jié)合多變的環(huán)境變量(溫度、降水等)進(jìn)行高海拔地區(qū)森林優(yōu)勢樹種識別需要進(jìn)一步研究。