胡圣元，龐 勇，蒙詩櫟，岳彩榮

(1. 西南林業(yè)大學(xué)林學(xué)院，云南昆明 650224；2. 中國林業(yè)科學(xué)研究院資源信息研究所，北京 100091；3. 國家林業(yè)和草原局林業(yè)遙感與信息技術(shù)重點實驗室，北京 100091)

森林變化包括森林突發(fā)性變化和時間序列趨勢變化[1]。森林?jǐn)_動是森林突發(fā)性變化，是研究者們在森林變化中關(guān)注的重點。隨著大量中高分辨率衛(wèi)星遙感數(shù)據(jù)的免費開放，研究者可利用美國Landsat系列、歐洲Sentinel系列和我國高分系列光學(xué)衛(wèi)星數(shù)據(jù)開展時間序列森林變化分析研究。目前基于時間序列的森林變化檢測方法大致包括3類：（1）基于多時期分類圖斑的方法[2]；（2）基于多時序像元級閾值判定的方法[3]；（3）基于長時間序列軌跡分析的方法。基于多時期分類圖斑的變化分析方法是通過對不同時期衛(wèi)星影像進(jìn)行地物識別和分類，并對分類結(jié)果進(jìn)行比較，從而確定森林變化的發(fā)生情況[2]。因此，時間序列數(shù)據(jù)質(zhì)量與分類精度都會對變化檢測的結(jié)果產(chǎn)生較大的影響?；诙鄷r序像元級閾值判定的變化分析方法是通過構(gòu)建植被指數(shù)，計算時間序列中逐個像元的指數(shù)曲線，在指數(shù)超出閾值時即可檢測到森林變化。這類森林變化分析方法包括植被變化追蹤算法（Vegetation Change Tracker, VCT）[4]、Hilker等[5]使用的映射反射率變化的時空自適應(yīng)算法，相比多時期分類圖斑的變化檢測，該類方法能夠在長時間范圍內(nèi)連續(xù)觀測同一地區(qū)的森林覆蓋情況，可以得到變化區(qū)內(nèi)更細(xì)致的時空變化特征，例如森林在某段時間內(nèi)經(jīng)歷的成熟-砍伐-恢復(fù)-成熟等活動。這類方法的不足是閾值的選取本身會影響變化檢測的結(jié)果，并且構(gòu)建的指數(shù)會受到輻射差異的影響[1]?；陂L時間序列軌跡分析的變化分析方法是通過時間序列光譜軌跡分析來檢測森林變化：Hansen等[6]使用的方法是曲線擬合法，即預(yù)先確定每種類型森林變化曲線的形狀，根據(jù)時間序列曲線形狀進(jìn)行變化歸類；Land-Trendr方法[7]是分割法，首先將時間序列光譜軌跡分割成一系列直線段以模擬時間軌跡重要特征，再根據(jù)分割所得直線段端點提取變化信息；CCDC方法[8]是統(tǒng)計邊界法，首先將時間序列曲線分解為趨勢、季節(jié)變化和噪聲，然后對時間序列的統(tǒng)計學(xué)邊界進(jìn)行估計，如果某一年的數(shù)值超過了估計的邊界值則認(rèn)為有森林變化發(fā)生。軌跡分析法相比閾值法，更好地利用了時間序列信息。

VCT算法是多時序森林變化檢測方法中的經(jīng)典方法之一，其構(gòu)建的IFZ指數(shù)（Integrated Forest Z-score，IFZ）簡便有效，運算速度較快。Huang等[4]使用VCT算法對美國東部7個州的森林動態(tài)變化進(jìn)行了檢測，總體精度達(dá)到80%。張連華等[9]通過纓帽變換構(gòu)建擾動指數(shù)（Disturbance Index）以代替IFZ進(jìn)行變化識別，但仍有檢測精度較低的年份。黃春波等[10]在IFZ的基礎(chǔ)上引入了NDVI進(jìn)行加權(quán)計算，對三峽大壩附近森林進(jìn)行變化檢測，相比Huang等在美國的結(jié)果，漏檢率和過檢率分別降低到了4.7%和5.3%。VCT算法對森林覆蓋變化劇烈的情況較為靈敏，如大規(guī)模砍伐、森林火災(zāi)等導(dǎo)致的森林?jǐn)_動等；對森林覆蓋變化輕微的情況如擇伐、少量病蟲害導(dǎo)致的森林?jǐn)_動以及森林逐漸恢復(fù)的檢測效果欠佳[4]。同時，有植被覆蓋的非森林區(qū)（如農(nóng)田、草地等）的變化過識別會影響森林區(qū)域的變化精度。對于農(nóng)田區(qū)，VCT算法利用擾動變化前后數(shù)年的信息以排除季相因素導(dǎo)致的偽變化[3]。

通常森林的變化情況在一段時間內(nèi)相對穩(wěn)定，較少出現(xiàn)短時間內(nèi)發(fā)生連續(xù)擾動，因此VCT算法能較準(zhǔn)確地檢測出森林發(fā)生擾動的年份。然而對于短輪伐期的速生用材林，森林砍伐后的再造林或萌生林相隔時間較短，僅經(jīng)過兩年甚至一年的恢復(fù)，IFZ指數(shù)便能恢復(fù)到正常的森林水平，則VCT算法易判定為森林沒有發(fā)生擾動，容易出現(xiàn)漏檢情況。因此，對IFZ指數(shù)進(jìn)行優(yōu)化以提高其算法對森林?jǐn)_動的檢測能力，可以減少偽變化同時降低漏檢風(fēng)險。

本研究以高峰林場為研究試驗區(qū)，在減弱VCT偽變化判讀機制的前提下使用了結(jié)合NDVI的IFZ指數(shù)進(jìn)行森林?jǐn)_動的檢測，對比使用原始IFZ指數(shù)的森林?jǐn)_動檢測結(jié)果，新的IFZ指數(shù)在保證生產(chǎn)者精度不變的前提下提升了用戶精度和總精度。

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

1.1 研究區(qū)概況

以廣西高峰林場銀嶺分場至延河分場為研究區(qū)，高峰林場位于南寧市興寧區(qū)北部，坐標(biāo)范圍22°50′～23°4′ N，108°8′～108°32′ E；屬濕潤的亞熱帶季風(fēng)氣候，陽光充足，雨量充沛，年平均氣溫在21.6 ℃左右，年均降水量達(dá)1 304.2 mm。森林經(jīng)營面積8.7萬hm2，森林蓄積量570萬m3，其中3.3萬hm2林地位于南寧市北部丘陵地帶，森林覆蓋率大于80%。主要樹種包括杉木(Cunninghamia lanceolata (Lamb.) Hook.)、馬尾松(Pinus massoniana Lamb.)、巨尾桉(Eucalyptus grandis×urophylla)、八角(Illicium verum Hook. f.)、毛竹(Phyllostachys edulis (Carr.) H. de Lehaie)、火力楠(Michelia macclurei Dandy)等[11]。高峰林場擁有大量短輪伐期的速生用材林，森林經(jīng)營活動頻繁，由砍伐產(chǎn)生的森林減少和恢復(fù)產(chǎn)生的森林增加較多，呈現(xiàn)出較頻繁的森林變化。

1.2 遙感數(shù)據(jù)獲取與數(shù)據(jù)處理

采用的遙感影像為Landsat8 OLI地表反射率產(chǎn)品，該產(chǎn)品已經(jīng)過大氣校正、正射校正、輻射定標(biāo)等處理，可通過美國地質(zhì)調(diào)查局地球資源觀測與科學(xué)中心官方網(wǎng)站下載（http://earthexplorer.usgs.gov）。行列號為125/44的影像完整覆蓋高峰林場研究區(qū)，研究選取了成像時間為2013年至2019年期間生長季（3—10月）的所有影像數(shù)據(jù)。根據(jù)云的光譜特性及衛(wèi)星影像的像元質(zhì)量評估信息等，使用影像合成法利用多期衛(wèi)星影像的非云像元，生成研究區(qū)2013—2019年逐年度的無云影像[12-13]。

1.3 地面數(shù)據(jù)獲取

外業(yè)調(diào)查于2018年5月11—23日進(jìn)行。對選取部分研究區(qū)的擾動圖斑樣點進(jìn)行實地調(diào)查，包括研究區(qū)邊界外部靠近研究區(qū)邊界的樣本在內(nèi)共采集到實測森林變化情況樣本點190個。對非森林區(qū)的樣本進(jìn)行系統(tǒng)采樣，在研究區(qū)及附近以橫向30 km、縱向21 km的矩形范圍進(jìn)行了間隔3 km的系統(tǒng)抽樣，選取了88個樣本點，利用衛(wèi)星影像對樣本點所處或距樣本點較近的擾動圖斑進(jìn)行目視判讀，獲取擾動信息。共采用278個樣本點對擾動檢測結(jié)果加以驗證。實地調(diào)查樣本點及系統(tǒng)抽樣目視判讀樣本點空間分布如圖1所示。

圖 1 外業(yè)調(diào)查樣本點空間分布圖Fig. 1 Distribution of field plots

2 基于綜合森林指數(shù)檢測森林?jǐn)_動

2.1 VCT算法介紹

VCT算法的基本原理是在影像中自動提取森林樣本，以此森林樣本作為參照，計算每個像元與純凈森林像元的相似性程度，即森林指數(shù)（FZ），再結(jié)合多波段或者多指數(shù)的FZ形成綜合森林指數(shù)（IFZ）。IFZ可作為每個像元是否為森林像元的判斷依據(jù)，根據(jù)每個像元IFZ指數(shù)的變化情況來判斷該像元在時間序列中所處的時間年份是否發(fā)生了森林?jǐn)_動[3]。

由森林像元的光譜曲線特征可以看出，其在紅光波段的反射率較低，當(dāng)森林樣本達(dá)到一定數(shù)量時，就會在紅光波段反射率直方圖的低值區(qū)域形成峰值，即森林峰值。通過提取該峰值區(qū)域的像元，同時對水體和陰影進(jìn)行掩膜處理，即可得到純凈的森林樣本像元[14]。對于整景Landsat影像，不同地區(qū)的森林可能由于物種類型、群落結(jié)構(gòu)、地形等因素導(dǎo)致的光譜差異，呈現(xiàn)略有差別的光學(xué)特征，若對全景影像采用直方圖森林峰值法提取森林樣本，會導(dǎo)致提取到的森林樣本較為單一[15]，因此可采用一定尺寸窗口將影像劃分為若干區(qū)塊，每個區(qū)塊單獨進(jìn)行直方圖森林峰值的提取。設(shè)置窗口時不宜過小，否則可能會因窗口未包含森林樣本而將農(nóng)田等其他地物的峰值作為樣本提取，從而導(dǎo)致森林樣本不準(zhǔn)確。在本研究中采用400×400像元窗口大小。

森林指數(shù)是基于像元統(tǒng)計特征的指數(shù)，它表示了每個像元與純凈森林像元的相似程度。FZ與IFZ的計算公式如下：

其中NB為參與運算的波段數(shù)目，bp為每一個像元在某一波段上的值，為i波段上森林樣本像元值的均值，SDi為i波段上森林樣本像元值的標(biāo)準(zhǔn)差。在使用單波段進(jìn)行計算時（NB=1），IFZ等同于森林指數(shù)FZ；使用多個波段時（NB＞1）即為綜合各波段的FZ值得到整合值IFZ。

由公式（1）可知，像元的IFZ值越低，其與純凈森林樣本越相似，越有可能是森林，反之則為非森林。排除水體和陰影對IFZ變化的影響后，每個像元IFZ值的變化情況可用于檢測發(fā)生森林?jǐn)_動的時間點。若像元點在某年的IFZ值與其前一年相比有顯著的提升，那么它就有可能是一個由森林向非森林轉(zhuǎn)變的變化（擾動），若像元點在某年的IFZ值與其前一年相比有顯著的降低，則可能是從非森林轉(zhuǎn)變到森林。因此設(shè)置衡量IFZ變化量的閾值d1作為判斷森林?jǐn)_動的依據(jù)，若前后年份IFZ的差值大于d1，則認(rèn)為極有可能發(fā)生擾動。但只有差值并不能確定發(fā)生了森林?jǐn)_動，因為非森林像元的IFZ波動同樣可以造成較大的IFZ差值。因此，設(shè)置一個閾值d2用于界定變化前后地物類型，只有當(dāng)潛在變化點的前一時刻的IFZ值低于d2時，這樣的變化才是一個由森林向非森林轉(zhuǎn)變的變化。

為了提高IFZ對森林與非森林植被的區(qū)分度，以減少檢測到的偽變化、降低過檢測率、提高用戶精度，本研究采用將NDVI與原始波段加權(quán)結(jié)合的方式，從而得到更好的森林?jǐn)_動檢測結(jié)果。研究流程如圖2所示。

圖 2 綜合森林指數(shù)檢測森林?jǐn)_動流程Fig. 2 Workflow of forest disturbance detection using forest z-scores

2.2 融合NDVI的IFZ森林指數(shù)

在已有的使用VCT算法的森林變化檢測方法中，參與運算IFZ指數(shù)的波段并不相同，Huang等人對Landsat 5 TM影像選用了紅波段和兩個短波紅外波段[3]；Pang等使用Landsat 7 ETM+影像時選取了藍(lán)、綠、紅、近紅外和兩個短波紅外共6個波段[16]。但有研究表明，雖然森林冠層在近紅外光譜區(qū)間的反射率很高，但是非森林的地表根據(jù)其覆蓋情況不同，在近紅外光譜區(qū)間的反射率或高或低，從而導(dǎo)致在這個波段上發(fā)生的反射率在某一特定方向上的變化并不一定代表森林覆蓋發(fā)生變化，因此不宜使用近紅外波段作為IFZ指數(shù)的計算因子[4]。此外，藍(lán)、綠、紅波段之間有著顯著的相關(guān)性[17]，因此也不宜同時使用這3個波段。植被指數(shù)也可以參與IFZ的計算。部分農(nóng)田的光譜與森林相近，其IFZ值也相對較低，有時可以和IFZ值較高的森林相當(dāng)，因此農(nóng)田等非森林植被會對森林?jǐn)_動檢測產(chǎn)生干擾，在高峰林場這樣變化頻繁的區(qū)域，部分森林?jǐn)_動產(chǎn)生的IFZ變化與農(nóng)田的IFZ變化更難以區(qū)分；將NDVI作為計算因子之一參與到IFZ的運算中，所得到的新的指數(shù)NIFZ，能較好地減弱農(nóng)田對森林?jǐn)_動檢測的影響[10]。

2.3 IFZ的波段加權(quán)

由IFZ的公式可知，像元的IFZ值是參與計算的每個波段FZ值的平均。在引入了NDVI作為新波段參與運算后，由于植被指數(shù)波段的取值范圍與反射率波段的取值范圍不同，例如可見光波段的反射率實際取值范圍大部分集中在0.2以下，NDVI值則較少出現(xiàn)低于-0.4的情況，因此不同波段的FZ存在取值范圍不等的問題。為了統(tǒng)一各波段FZ值的取值范圍，需將各波段FZ結(jié)果按比例縮放到同一尺度。將參與計算的波段乘以同一個放大系數(shù)，從而將各個波段FZ換算成具有相同權(quán)重的指數(shù)因子。因此，本研究將參與運算的首個波段作為基準(zhǔn)，以各波段FZ值的均值作為確定縮放比例的參考，即為NIFZ：

2.4 計算流程與驗證方法

使用合成影像的藍(lán)、紅、近紅外和兩個短波紅外波段進(jìn)行計算和分析。首先使用NDVI與近紅外（845～885 nm）、藍(lán)波段（450～515 nm）對水體和陰影進(jìn)行閾值法掩膜；接下來對影像的紅波段（630～680 nm）以400×400像元大小為窗口遍歷影像，利用直方圖的森林峰值提取森林樣本，得到每一影像的森林樣本，并以此樣本計算NIFZ，以及由此衍生的NIFZ2和NIFZ3。最后，計算紅、短波紅外1（1 560～1 660 nm）和短波紅外2（2 100～2 300 nm）3波段的IFZ指數(shù)進(jìn)行對照。

為了對比各指數(shù)在區(qū)分森林像元和非森林植被像元的表現(xiàn)，分別選取了不同年份、不同位置的若干森林區(qū)和非森林（農(nóng)田為主）像元的4種指數(shù)值進(jìn)行了對比。從圖3可以看出森林像元的NIFZ2與NIFZ3值的分布比IFZ和NIFZ更為集中，而就逐年取值分布來看，NIFZ與NIFZ2的森林與非森林像元的區(qū)分度更大。4種指數(shù)在森林與非森林之間都存在著重疊區(qū)間，相比之下NIFZ2的重疊區(qū)間更小。從曲線形態(tài)來看（圖4），4種指數(shù)的曲線形狀是相近的，都能夠反映森林?jǐn)_動變化特征。此外，根據(jù)4種指數(shù)逐年取值分布情況可以選定森林像元與非森林像元指數(shù)值的區(qū)分閾值d2，進(jìn)行后續(xù)的森林?jǐn)_動檢測的運算。在此，對IFZ、NIFZ、NIFZ2和NIFZ3，d2的值分別取4.5、3.5、2.5和2.3。

在分別利用各指數(shù)得到對應(yīng)的森林?jǐn)_動檢測結(jié)果后，利用驗證點擾動信息與對應(yīng)位置的檢測信息進(jìn)行對比，逐年記錄檢測情況與實際情況，統(tǒng)計各指數(shù)檢測結(jié)果的漏檢測率、過檢測率和總精度，以進(jìn)行各指數(shù)的精度評價。

圖 3 4種綜合森林指數(shù)在森林和農(nóng)田的取值區(qū)間對比，末尾的f與n分別表示森林和農(nóng)作物Fig. 3 Comparison of value ranges of 4 z-scores in forest(-f) and crop(-n) areas

3 結(jié)果與討論

采用IFZ、NIFZ、NIFZ2和NIFZ3作為反映森林?jǐn)_動指數(shù)應(yīng)用于VCT算法中，在不使用前后年份進(jìn)行驗證的情況下得出森林?jǐn)_動檢測的結(jié)果，并使用外業(yè)調(diào)查數(shù)據(jù)結(jié)合衛(wèi)星影像來進(jìn)行二者的精度檢驗。其中NIFZ2得到的結(jié)果如圖5所示。由不同指數(shù)擾動檢測結(jié)果對比圖（圖6）來看，4個指數(shù)的表現(xiàn)較為相似，但NIFZ與NIFZ2的檢測結(jié)果相比IFZ的結(jié)果減少了很多細(xì)碎的圖斑?？傮w而言，NIFZ2的檢測結(jié)果減少了大量主要類型為2015年森林?jǐn)_動的孤立像元，從檢測結(jié)果圖上看是效果最好的指數(shù)。

最終，使用采集和抽樣的驗證點對4種指數(shù)的檢測結(jié)果進(jìn)行驗證，各指數(shù)的檢測精度如表1所示?？傮w情況顯示，NIFZ2的精度最優(yōu)，總精度達(dá)79.2%，NIFZ其次，為77.4%；而從年度情況來看，2014年4個指數(shù)的用戶精度都比較低，而2015年NIFZ2的用戶精度獲得了較大提升，其余3個指數(shù)的用戶精度雖比2014年有所提升但仍然不理想。精度較低主要是由于2013年和2014年的影像中殘存的云像元對檢測結(jié)果有一定影響；3種指數(shù)中NIFZ3的精度最低，說明對不同年份的森林指數(shù)結(jié)果進(jìn)行加權(quán)比較并不能消除不同年份之間的差異，同時可能還會對森林?jǐn)_動的識別帶來負(fù)面的影響，這可能是由于盡管不同年份之間森林像元指數(shù)值存在差異，但利用森林峰值提取森林樣本時，不同年份的影像提取到的森林峰值的反射率值較為接近，所以被選為森林樣本的像元的指數(shù)值的差異在不同年份之間并不大，因此對不同年份加權(quán)并沒有令NIFZ3在相同地物上不同年份間的表現(xiàn)趨近相同。

圖 4 4種綜合森林指數(shù)在擾動區(qū)的時間序列特征Fig. 4 Curves of the four z-scores in a forest disturbance area

圖 5 應(yīng)用NIFZ2的VCT森林?jǐn)_動檢測Fig. 5 VCT forest disturbance detection using NIFZ2

圖 6 應(yīng)用4種綜合森林指數(shù)的VCT森林?jǐn)_動檢測局部對比，對應(yīng)圖5白框Fig. 6 Comparison of VCT forest disturbance detection results of four z-scores, located in the white box of Fig. 5

表 1 4種綜合森林指數(shù)變化檢測的精度對比Table 1 Comparisons of change detection accuracy using four different z-scores

高峰林場的驗證結(jié)果表明，采用NDVI的3個指數(shù)所得到的森林?jǐn)_動變化檢測結(jié)果的精度均好于采用原始影像光學(xué)波段的IFZ的擾動檢測結(jié)果。相比IFZ的結(jié)果，3種使用NDVI并加權(quán)的指數(shù)的用戶精度得到了顯著提升，而生產(chǎn)者精度則有小幅度的降低。由于Landsat8地表反射率產(chǎn)品未經(jīng)過地形校正，處于向陽坡面的森林，其指數(shù)取值會比背陽面更高，更加接近農(nóng)作物等非森林植被的指數(shù)；此外，不同質(zhì)的森林由于樹種、密度、郁閉度等差異，也會影響森林指數(shù)值。這些都是利用閾值判斷的森林變化檢測方法的不確定因素，加權(quán)融合NDVI與IFZ的方法增加了森林與非森林在森林指數(shù)上的區(qū)分度，間接降低了這些不確定因素的影響，從而提高了森林?jǐn)_動檢測的精度。其他遙感衍生的植被指數(shù)例如增強型植被指數(shù)（EVI）和差值植被指數(shù)（DVI）也能夠反映森林覆蓋變化相關(guān)的信息，但對于上述指數(shù)是否適用于提高IFZ的檢測精度有待進(jìn)一步驗證；此外，經(jīng)過穗帽變換所得到的波段，同樣含有森林變化的信息[18]，也有著與IFZ相結(jié)合的可能性。

4 結(jié)論

通過將植被指數(shù)引入時間序列的綜合森林指數(shù)，提供了植被變化追蹤算法在森林?jǐn)_動頻繁區(qū)域的適應(yīng)性，主要結(jié)論如下：（1）在NDVI與IFZ的結(jié)合方式中，采用森林樣本FZ均值加權(quán)的NIFZ2指數(shù)的檢測精度最優(yōu)，其精度比IFZ提升了12%；（2）使用年度森林樣本FZ均值對不同年份指數(shù)進(jìn)行加權(quán)并不能有效提高擾動檢測的精度；（3）通過在IFZ中加入NDVI因子，使算法的誤檢率大大降低，平衡了快速變化區(qū)域內(nèi)的漏檢與誤檢的問題。